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6.1 



ABSTRACT 



A least squares algorithm is developed to solve for the 
trajectory and transponder array coordinates of the current 
velocity profiler, Pegasus. Measurement residuals and 
parameter precision are computed for data quality analysis. 
Travel times from a maximum of four seafloor transponders, 
pressure sensor depths, and transponder positions are input 
with their respective accuracy estimates. The algorithm is 
used to analyze a 2250 m profile from the Monterey Canyon with 
four transponders, one of which had not been positioned. This 
transponder's unknown position is found and problems in the 
other array coordinates identified. Transponder coordinate 
precision improves by factors of ten in the horizontal and 
five in depth, to about 13 m (Drms) and 2 m (lo) respectively. 
Trajectory precision is about 7 m (Drms) horizontal, with high 
correlation between points. Thus, the precision of horizontal 
velocity components, determined by time differencing points, 
is better than 11 cm/s (lo) . Depth precision is better than 
3m (lo), except in the deeper portions where anomalous 
pressure residuals near the depth of the transponder array 
suggest systematic pressure errors needing further study. 
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I. 



INTRODUCTION 



The Oceanography Department of the Naval Postgraduate 
School is presently conducting a study of ocean current 
velocities in the waters off the central California coast in 
the area of the Monterey Canyon. During the course of the 
study, several stations have been established for the 
collection of data. Station 10, northernmost in this chain of 
stations and the primary focus for analyses presented in this 
thesis, contains four seafloor transponders, one of which was 
considered to have an essentially unknown position at the 
beginning of the work described here. 

The main oceanographic instrument used in this survey is 
the dropsonde Pegasus. As it descends and ascends through the 
water column, it records at 16 second intervals, both the data 
from oceanographic sensors and the time for return travel of 
acoustic signals from transponders previously set in place on 
the ocean floor. 

In the past, Pegasus-generated velocity studies have used 
the raw data in various ways. The trajectory of the 
instrument has been fitted to a curve which was then differen- 
tiated with respect to time in order to produce the horizontal 
velocity components (Halkin et al., 1985). Alternatively, and 
as is the present practice at NPS, velocities have been 
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calculated from position differences and then smoothed with a 
spatial filter such as a seven point running average. 

This study outlines a method of determining the current 
velocities and their precision by applying a least squares 
adjustment procedure, adapted from geodetic survey techniques, 
whereby all observational data are used simultaneously. 
Specifically, three basic problems will be analyzed: 

* The precision of the transponder coordinates from which 
the position of Pegasus is derived. 

* The precision of the Pegasus positions and consequently 
the current velocities derived from Pegasus navigation. 

* The depth at which pressure becomes critical in solving 
for velocities. 

Discussions of this topic will appear in the following 
sequence : 

* Chapter II provides background information on the 
transponder network and the Pegasus system and its 
positioning . 

* Chapter III discusses sources of error resulting from 
systematic and observational inaccuracies in the data. 

* Chapter IV gives a brief explanation of the least squares 
adjustment process. 

* Chapter V describes the final procedure used. 

* Chapter VI discusses the results obtained. 

* Chapter VII provides conclusions and recommendations for 
future consideration. 

* Appendices contain algorithms and explanations of the 
Fortran programs used. 
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II . BACKGROUND 



The area being surveyed is composed of ten stations 
located just off the California coast approximately between 
36° 05' and 36° 39' North Latitude, and 122° 08' and 124° 13' West 
Longitude. Of particular interest here is Station CIO, 
northernmost in the chain and located in the Monterey Canyon 
about ten miles west of Monterey (see Figure 1) . 

A. TRANSPONDER NETWORK 

At this station four transponders were deployed, horizon- 
tally separated by distances somewhat equivalent to transpon- 
der depths (around 2000 m) . Three of these transponders, 
responding at 11.5 kHz, 12.0 kHz, and 12.5 kHz, respectively, 
were located by the survey vessel according to traditional 
methods as described below. However, the bandwidth of the 
shipboard receiving instrument was too narrow to pick up the 
13.5 kHz signal from the fourth transponder, and therefore its 
position was largely unknown. (Table 1 defines transponder 
abbreviations, their frequencies, and their approximate 
depths. Figure 2 shows the survey net for Station CIO.) 
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Figure 1. General Area of Station CIO 
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TABLE 1 



CIO TRANSPONDERS 



Transponder 


1 


T1 


12 . 5 


kHz 


1880 


in 


Transponder 


2 


T2 


12.0 


kHz 


1990 


m 


Transponder 


3 


T3 


11.5 


kHz 


2230 


m 


Transponder 


4 


T4 


13.5 


kHz 


1890 


m 




Figure 2. Station CIO Survey Net 
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Depth (m) 



B. SURVEYING A TRANSPONDER NETWORK 



1 . Transponder Depth 

To determine each available transponder depth, the 
survey ship homes in on the instrument's signal until a 
minimum time is recorded, thereby assuming the ship to be 
directly over the instrument (see Figure 3). An appropriate 
harmonic mean sound velocity multiplied by one way signal 
travel time produces the depth of the transponder. 



SH I P 
TRACK 



— > 




Figure 3. Determining Transponder Depth 
2 . Baselines 

To establish the relative horizontal x and y coordi- 
nates of transponders, baselines are determined between the 
instruments. Running at a relatively slow speed (say 4-6 
kts) , the survey vessel attempts to cross each baseline at a 
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90“ angle. The process is repeated in the opposite direction. 
This procedure is then duplicated, several more times if 
possible. Careful observation of the analog trace that 
records round-trip travel times of the acoustic signal between 
the survey vessel and the two transponders indicates when the 
shortest distance has been reached, as shown in Figure 4. 



BASCLINE CRDSSING 




Figure 4. Analog Trace of a Baseline Crossing 

Signals show up on the trace as parabolas approaching 
one another. They will lie vertically in line if the baseline 
is crossed at a 90“ angle and therefore the shortest distance to 
each of them is reached simultaneously. At this point the 
ship will lie in a vertical plane containing the two 
transponders. The sum of the two horizontal ranges will be 
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at a itiinimuin and therefore establish the baseline (see Figure 
5) . 




Figure 5. Baseline Determination 

The time of crossing, ship's course, latitude/ 
longitude, and Loran coordinates are noted. Time units for 
each minimum range are measured and then multiplied by the 
appropriate harmonic mean sound velocity to obtain slope 
distances. These distances, the depths of the two transpon- 
ders, and the ship's course are then used to calculate the 
horizontal baseline length and the azimuth between the two 
transponders. Finally, a local xy plane coordinate system is 
established, setting the position of one of the transponders 
at 0,0, and positioning the other transponder relative to it. 

Many factors combine to make the information on this 
analog trace subjective and inaccurate. The analog paper may 
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be set to progress at different rates producing sweeps of 
0.5 s, 1 s, or 2 s, etc. The trace may actually "wrap around" 
one or more times as the time interval gets longer and longer, 
making page integer resolution difficult. At a one second 
sweep, the full page represents (after conversion from time to 
distance) meters traveled in one second, i.e., 1500 m/s. 

Furthermore, if the baseline is not crossed perpendic- 
ularly, a time offset between minimum ranges is observed which 
decreases the length of the baseline. This error can be 
removed to some extent by multiplying the time offset by the 
ship's speed and then by geometric relationships, computing 
the correction. This situation is illustrated in Figure 6 
where T1-T2 represents the actual baseline. A, B represent 
respective points on the analog trace where minimum slant 
ranges are indicated, and a, b represent respective minimum 
slant ranges to each transponder. 




SHIP TRACK 




Figure 6. Time Offset 
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In addition to the above possibilities for error, an 
inaccurate course heading at the time of baseline crossing 
will alter the geometry of the array, resulting in inaccurate 
coordinates. Further discussions of transponder depth and 
baseline determinations can be found in Kuo (1985), McKeown 
(1975), and Hart (1967). For an assessment of transponder 
network accuracy, refer to Chapter III. 

Figure 7 illustrates the CIO transponder network and 
the approximate position of Pegasus drop #157 positioned in a 
local coordinate system using T1 (12.5 kHz) as 0,0. 







Figure 7. The CIO Local Coordinate System 
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C. DESCRIPTION OF THE PEGASUS SYSTEM 

The oceanographic instrument Pegasus is a free-falling, 
acoustically-tracked velocity profiler. It was developed by 
H.T. Rossby and D. Dorson at the University of Rhode Island a 
decade ago to fill a need for an instrument which would accu- 
rately record the fine scale vertical structure in the ocean 
and prove inexpensive and easy to handle at sea. It provides 
a means to measure absolute velocity components throughout the 
water column. (For the seminal paper on Pegasus, see Spain et 
al., 1981.) Prior uses of Pegasus in studies of the Gulf 
Stream off the east coast of the U.S. have been documented in 
a paper and data report by Halkin et al . (1985). 

The position of the instrument as it descends and ascends 
is determined by travel time of signals emitted by Pegasus at 
10 kHz every 16 seconds which are heard and then answered at 
other frequencies by transponders previously set in place on 
the ocean floor, and by pressure readings. These time 
intervals, and temperature, conductivity, and pressure data 
from oceanographic sensors, are recorded and stored in a 
microprocessor-controlled memory in the instrument to be later 
down-loaded aboard ship when the Pegasus is recovered. 

The Pegasus model currently adapted for use by the Naval 
Postgraduate School is manufactured by Benthos. It is housed 
in a 17-inch glass sphere protected by a hard hat and, 
according to manufacturer claims (Benthos, 1989), provides 
operation to full ocean depth, integral flotation, no 
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corrosion, and large battery capacity allowing for 100 deploy- 
ments without opening. (A normal NPS school cruise uses about 
20 deployments.) It weighs approximately 45 kg in air and 
carries eight kg of expendable weight which sink the profiler 
in this present study at about 27 m/min. Six receiver 
channels for communication with transponders are available 
with frequencies located at 0.5 kHz intervals ranging from 
11.5 kHz to 14.0 kHz. NPS supplied SEA BIRD model SBE-3 
temperature sensor and SEA BIRD model SBE-4 conductivity 
sensor are externally mounted with electrical interfaces, 
along with a Paroscientif ic model 410KT pressure 
transducer with a companion Intelligent Transmitter circuit 
board . 

D. PEGASUS POSITIONING 

1 . Computation of Pegasus Coordinates 

The position of Pegasus at each 16 second interval can 
be determined by solving a system of equations using the 
formula: 



Time = Distance/Sound Velocity 

Using one-way signal times from two transponders (see 
Figure 8) , the depth of Pegasus as derived from pressure 
information, and an appropriate sound velocity, there will 
then be two equations in two unknowns, allowing for a solution 
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PEGASUS 




Figure 8. Determining Pegasus Position 



of XP and YP at a given Pegasus position. The solution is 
derived from the equations: 



Time^ 



Time2 



[ (XP - XJ^ + (YP - YJ^ + (ZP - Zi)^]^'Vv 
[ (XP - Xa)^ + (YP - Y2)^ + (ZP - Z2)^]'^VV 



where : 



Time^ 



one-way travel time from transponders 1,2; 



V = effective sound velocity along the path between 

Pegasus and the transponder; 

XP,YP,ZP = Pegasus position coordinates (ZP is assumed 
known from the simultaneous pressure 
observations) ; 
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Xi,Yi,Zi = transponder coordinates (assumed known). 

The traditional oceanographic method outlined above 
calculates Pegasus position coordinates. However, it lacks a 
means of obtaining precision estimates for the derived 
positions, and it is not able to use all the observed data 
simultaneously. Furthermore, there is no redundancy in the 
computational process and thus systematic errors or blunders 
can escape unnoticed. 

2 . Estimating Horizontal Current Velocity 

Current velocity vectors u and v (describing velocity 
in the x and y directions, respectively) between two Pegasus 
positions, can be obtained by taking the difference between 
coordinates and dividing by the time interval: 

u = (XP2 - XPi)/dt 

V = (YP2 - YPi)/dt 

dt = 16 s 

This method, like the one mentioned in the previous 
section, does not provide estimates of precision and does not 
make use of all available data. 
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III. SOURCES OF ERROR 



Davis et al . (1981, pp. 15-20) defines measurement error 
as the difference between a measured and a "true" value and 
describes these errors as being generally of three types: 

* blunders or mistakes. 

* systematic errors. 

* random errors. 

Blunders can be caused by carelessness, equipment failure, 
or false interpretation. They are usually large enough to be 
easily spotted when results are analyzed. 

Systematic errors follow a defined pattern or system, 
which when discovered, can usually be expressed 
mathematically. Such factors as observer limitations, 
instrument imperfections or inadequate calibration, 
meteorological conditions, or poor choice of mathematical 
model can produce the pattern which will remain consistent as 
long as the elements of the system remain the same. 
Systematic errors are not eliminated by repetition of 
measurements. They must be ferreted out and either removed 
from the observations or their effects added to the model. 
Thus, to begin such an analysis, serious consideration has to 
be given to locating and defining errors resulting from 
systematic and observational inaccuracies in the data. 
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Random errors are those that remain after mistakes and 
systematic errors have been accounted for. The values of 
these errors should be unbiased and will ideally distribute 
themselves about a mean of zero. It is these random errors 
that the least squares process attempts to minimize. In the 
analysis described in this thesis, a considerable amount of 
time was spent in determining what the estimates of the a priori 

standard deviations of these random errors should be, both 
from the point of view of reasonable physical reality and 
meaningful results. The adjustment technique uses these a priori 

estimates to weight the contribution of each observation to 
the final solution. 

The following potential error sources were studied and, 
where applicable, appropriate standard deviations of the 
errors were then entered as weights into the least squares 
adjustment model: 

* Transponder coordinates. 

* Signal time. 

* Velocity of sound. 

* Pressure/Depth relationship. 

A. TRANSPONDER SURVEY 

Of prime importance is the precision of the transponder 
positions. Determining these coordinates is a difficult and 
time-consuming process. As noted before, in the CIO 
transponder net, the location of the 4th transponder could not 
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be determined by traditional methods because the band width of 
the ship's receiver was not wide enough to pick up the 



13.5 kHz signal. Also, it was not possible at the time to 
determine more than two of the three baselines between the 
three "known" points, thereby precluding a mathematical 
closure of the triangle. 

1 . Transponder Depth 

To establish a reasonable estimate for the standard 
deviation of transponder depth error, it was noted that any 
horizontal offset resulting from the ship's not being directly 
overhead always produces a positive depth error, i.e., the 
slant range must be longer than the vertical. Figure 9 is 
similar to Spain et al. (1981, p. 1557) and illustrates this 
process . 



The solution is as follows: 



(H + h)2 = 



2Hh + h^ = 



2H 





where : 



H = true transponder depth 



h = depth error 
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Figure 9 . Transponder Depth Error 



If h << H, the last term can be disregarded, and 



h ~ 




A realistic estimate of the ability of a survey vessel 
to cruise directly over a transponder by using this method is 
about 100 meters (Schnebele, 1990a) . Given the uneven canyon 
terrain of this survey which makes resolving the position of 
the shoalest depth even more difficult, it seemed reasonable 
to ascribe an error of 200 m to the offset. Thus, if the true 
depth lies at 2000 m and there is a horizontal offset of 
200 m, there will be a depth error of about 10 m. 
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Therefore, a = 10 meters was chosen for transpon- 
ders 1, 2, and 3. Since there was no observed depth for 
Transponder 4, a depth of 1900 meters was read from the latest 
NOAA chart using the position of transponder drop^ and, since 
it had not been surveyed in, a larger standard deviation of 
100 meters attached to it.^ 



^ (T1z,T22,T32) ~ ^ 



= 100 m 



2 . Baselines 

In determining baselines, three sources of error 
predominate: 

* Minimum range reading taken when the ship is offset from 
the actual baseline. 

* Transponder depth error propagating into the computed 
baseline length. 

* Error in ship's course affecting baseline azimuth. 

In a simplified ideal situation where the baseline is 
crossed at a 90° angle in the center of the line and the 
depths of the two transponders are the same, the expected 



^The transponder position was calculated from Loran 
signals which are based on North American Datum 27. The NOAA 
chart uses NAD83. Positions may vary up to 100 meters between 
the two datums in this locality. 

^After the October 17, 1989 earthquake, a depth measure- 
ment of 1931 meters was made of the 4th transponder by NPS 
scientist Tarry Rago in a submersible. 
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error in baseline length from measured slant ranges can be 
calculated as follows (see Figures 10 and 11). 
a. Baseline Error Due to R 




N B H 



Figure 10. Baseline Error 

R = distance off baseline when slant ranges were 
measured; 

B = true baseline between transponders (B = + B 2 ) ; 

b = baseline error (b = bi + h^) ! 
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(Bi + bi)^ = 



2Bibi + bi^ = 



bi(2Bi + bi) = R2 



b, = 



R 



R 



2B^ +b^ 






R 

2B, 



(1 - 



2B 



1 2 

+ (^) 



■2B, 



Since + . . . ) « 1, 



R 



bi ^ and, similarly, bz ^ 



R 



- 2B 



- 2B, 



If, as stated above, B^ = Bz = B/2, then b^ = bz = b/2 



+b2 



2 2 
R ^ R 






2B^ 2B_ - 2 'B/2 B/2' - 



^(i) 

2 ^b' 



b„ ~ 



2R^ 



R - B 



This error increases the length of the baseline. 
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b. 



Baseline Error Due to h 



Bl > ^ 



D 




Figure 11. Relationship of Baseline and Depth Errors 

H = true transponder depth; 

h = depth error; 

S. = slant range from ship to transponders 1,2. 

+ Bi^ (1) 

= (Hi + hi)2 + (Bl + bi)2 (2) 

Si remains constant. 

Eq. (2) - Eq. (1) : 
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0 



= 2Hihi + + 2Bibi + 



2B^bj^(l +23^) +_^j 









But: 



(1 " 2 H 7 ) 



(1 " 217 ) 



^1^1 ^1 ^1^ 

p— = Bj^b^d - j— + (^— ) - ...) 

<1 "2ii7' 

-H^h^ b b^ 2 

b^= -Vl^l - 2B7 - •••) 






h. 



Since (^+ ...) « 1, and (^+ ...) « 1, then 

z H ^ 



Bjbj ^ “Hihi 



bi 



«ihi 

B, 



If Hi = H 2 = H, hi = h2 = h, Bi = B 2 = B/2, bi = b2 = b/2, 



b. = 



bi + b2 ^ -Hh(^ + 

®1 ®2 



-">'<5^ * 5^) 



bh ~ - 



4Hh 

B 



then 
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As shown in Figure 11, the hypotenuse representing 
the given slant range of the signal must be the same in both 
triangles ADE and ACF. FC must equal ED. Therefore, if the 
assumed depth is increased by error, then the computed 
baseline must consequently decrease, 
c. Total Baseline Error 

Thus the total baseline error b^ + b^ becomes: 



t>Total ~ 






4Hh) 



Using typical figures of: B = 2000 m, H = 2000 m, 
h = 10 m, and R = 100 m, it should be possible to obtain a 
value for bjotai = ±30 m. 

If the transponder depths are not the same, and 
the baseline is not crossed right at the center but is crossed 
on a perpendicular heading, the baseline error will be: 



-’Total ~ 






2 'B, 



B, 



H,h, H„h 

^ B, B., ^ 



In addition to baseline errors discussed above, an 
incorrect azimuth will add even more uncertainty. An azimuth 
error of 0.5“ over a baseline of 2.0 km, for example, will 
produce an error in position of one end of the baseline with 
respect to the other of 18 meters. 
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Based on the above analysis, it was initially 
assumed that the relative horizontal coordinates of Tl, T2 , 
and T3 would be known to approximately ±30 m. Hart (1967) , 
for example, claims that with careful repetitive procedures 
involving all baselines in a transponder array, relative 
accuracies of the order of five meters can be achieved. 

As preliminary data analysis proceeded with the 
CIO network, it became obvious that there were major inconsis- 
tencies with the given baseline data. Preliminary adjustments 
(see Chapter V) would not converge when the given transponder 
coordinates were used. Further analysis suggested both length 
and azimuth problems with the given T2-T3 baseline. The fact 
that the T3-T1 baseline had never been determined led to a 
situation in which no positional redundancy existed in the Tl, 
T2 , and T3 network. As a result of these problems, initial 
positions of Tl, T3 , and T4 were determined from the Loran 
coordinates of their drop sites, while T2 was computed 
relative to Tl from baseline survey data. 

As a conseguence of this less than ideal 
procedure, it was considered that standard, deviations of ±100 
m should be allocated to the transponder coordinates of T2 , 
T3 , and T4, while Tl would be assumed fixed in order to 
provide a positional datum (albeit a rather arbitrary one) for 
the array. The 100 m allocated for positional standard 
deviations seemed reasonable in the light of both the known 
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positional accuracy of Loran and the added uncertainty from 
drift as the transponders settled to the bottom. 

In summary, therefore, the following standard 
deviations were assumed: 



(T2 ,T3 ,T4 

x,y' x,y' X, 



100 m 



= 0.01 m 

x,y 



B. TIMING 

The time observations required a consideration of how well 
the acoustic time intervals could be established between 
Pegasus and the transponders. The travel times were assumed 
to be independent of one another. The Pegasus manual states 
that the transponders have a 12.5 ms output pulse delay 
accurate to within ±0.5 ms (Oceanographic Instrument Systems). 
For the initial transponder adjustment, = 0.002 s was 
chosen as a reasonable weight (see Chapter V) , and then later 
reduced to 0.0005 s for the full run with improved transponder 
coordinates . 



Ojimel = 0.002 s 



OTime2 = 0.0005 S 
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C. VELOCITY OF SOUND 

Fofonoff (1963) states that velocity of sound is a 
function of the thermodynamical and chemical state of the 
water and is determined by using any complete set of variables 
of state such as temperature, pressure, and salinity. He goes 
on to observe that the precision of sound velocity must 
consider not only that of the oceanographic data, but also of 
the empirical formula used to convert these measurements to 
sound speed. 

Pegasus provides temperature and conductivity data which 
is converted to a sound velocity profile for the full range of 
the drop by algorithms published in the Unesco Technical 
Papers in Marine Science #44 (Fofonoff and Millard, 1983, pp. 
11-12, 49). 

In dealing with slant ranges, as in this case, acoustic 
refraction should be considered as well. Spain et al. (1981, 
pp. 1562-1564) presents an equation derived by Vaas (1964) 
which corrects for the effects of ray bending, taking into 
consideration the relative positions of the projector and the 
receiver even when they are at similar depths. This equation 
is stated to have an accuracy of better than 0.2 m/s for all 
angles and depths. 

Another excellent source for acoustic refraction informa- 
tion in a survey situation is the SASS Accuracy Study 
Simplified Ray Bending Correction (General Instrument Corp. , 
1975) and its follow-up Ray Bending Correction for Depth 
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Sounders, An Informal Approach (General Instrument Corp. , 
1976) . These two papers develop the equations for an 
effective sound velocity in detail and include graphs showing 
the errors relative to lateral angles. 

In this study it was found that the use of an average 
harmonic mean sound velocity produced acceptable results. To 
obtain this average, the sound velocity profile was used to 
calculate harmonic sound velocity profiles with respect to 
each transponder. These four harmonic mean profiles were then 
averaged (see Figure 12) and the resulting profile used to 
calculate the approximate X and Y coordinate of each Pegasus 
position. The harmonic mean profiles differ for each 
transponder because of the wide difference in deployment 
depths. Harmonic means at the depths found here differed by 
about 1.0 m/s from the average used, representing a bias of 
less than 0.7 m in a typical one km path length. 

Refraction was not regarded as significant for this study 
because of the relatively short path lengths and small sound 
velocity gradients encountered (Schnebele, 1990b) . In deeper 
waters with longer paths, refraction may be significant. 

The harmonic mean used in this model, therefore, assumes 
no refraction. The small error that it introduces was felt to 
be inconsequential compared to the uncertainties in 
transponder coordinates and signal travel time measurements. 
Thus the harmonic means used in the adjustment were assumed to 
be exact quantities without significant random error. 
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Figure 12. Harmonic Sound Velocity 

D. PRESSURE/DEPTH RELATIONSHIP 

The algorithm used to convert pressure to depth is the 
Saunders and Fofonoff formula which uses the hydrostatic 
equation and the Knudsen-Ekman equation of state (Fofonoff and 
Millard, 1983, p. 25). The formula includes variation of 
gravity with latitude and depth and assumes standard sea 
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water. This formula is said to deviate by only 0.08 meters at 
5000 decibars from estimates based on EOS80, a considerably 
smaller error than those in pressure measurements as shown 
below. 

The manufacturer of the Pegasus pressure transducer claims 
a typical accuracy under difficult environmental conditions of 
0.02% (Paroscientif ic, 1986). At our full scale of 2200 m, 
this would compute to about 0.4 m. As the analysis 
progressed, further inconsistencies in the data began to lead 
to a suspected bias in pressure measurements, perhaps caused 
by a temperature hysteresis as Pegasus drops and rises through 
the deepest part of the water column. Therefore, a o =2.00 
m was taken to be the standard deviation of the Pegasus depth 
measurements, i.e.. 



o = 2.00. 

P„ 
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IV. LEAST SQUARES ADJUSTMENT 



The specific task of this study was essentially to analyze 
three basic problems: 

* The precision of the transponder coordinates, especially 
those of T4 , from which the position of Pegasus is 
derived. 

* The precision of the Pegasus positions and consequently 
the current velocities derived from Pegasus navigation. 

* The depth at which pressure becomes critical in solving 
for velocities. 

To provide a means of answering these questions, a Least 
Squares Adjustment Program was developed by Dr. John Hannah.^ 
This adjustment program, as well as a brief explanation, are 
provided in Appendix A. 

Least squares adjustment techniques are used to determine 
unique solutions for unknown parameters when there are 
redundant observations, i.e., more than necessary to specify 
the model (Davis et al . , 1981, pp. 38-39). The least squares 
estimator is an unbiased minimum variance estimator which is 
unique, mathematically easy to derive, and, when compared to 
other estimation techniques, leads to a smaller dispersion of 
random errors. It also is distribution free in the sense 



^Adjunct Research Professor, CNOC Chair in Mapping, 
Charting, and Geodesy, U.S. Naval Postgraduate School, 
Monterey, CA. , 1988-1990. 
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that a distribution is needed only for confidence interval 
testing (Hannah, 1989) . 

The mathematical model for the Pegasus data adjustment has 
the general form: 



Lai = Fi(xJ 



La2 ~ ^2 ( ) 



in which the first set of observations comes from Pegasus 
itself in the form of return signal travel times from the 
transponders to the instrument. The second set comes from an 
a priori knowledge of any of the system parameters such as the 

positions of the transponder coordinates, and each Pegasus 
depth as calculated from pressure. In general terms, any set 
of adjusted observations expressed as a function of a set of 
adjusted parameters can be written in a matrix expression: 

La = F(xJ 

The following is adapted from a least squares adjustment 
procedure written by Hannah (1981) and is used here with the 
author's permission. 

. . .These functions may be linearized by taking a first order 
Taylor series expansion about some approximate values for 
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the parameters, X 



The first function then becomes 



DF 

= TilT" [. , (>^a-Xo) + F(Xo) 



a X =X„ 
a 0 



or 



Lb> + Vi = AiX + Lo* 



or 



Vi = AjX + Li 



in which 



nr 



Ai = 






-I , X = (X,-Xo) 



and 



Li = Lo^ - 



Similarly, the second function may be linearized to give 

A 

V2 = A2X + L2 



In the above, Vj and the vectors of residuals arising 

from the observations Lb^ and Lb^ with the adjusted values of 
the parameters being given by the vector X^. Assuming that 
the two sets of observations are uncorrelated, then the 
least squares minimum variance estimate of X based on these 
two sets of observations is given by minimizing the function 

(}) = Vi^PiVi + V2^P2V2 - 2Ki^(Vi-AiX-Li) - 2K2‘^ ( V 2 -A 2 X-L 2 ) 

A 

With respect to the unknowns Vj, V2, Kj, K2, and X. The 
weight matrices for each set of observations are given by Pj 
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and ?2, in which Pi = E,“i and Pj = E2 “ 1 / i.e. the inverse of 

^ b ^ b 

the var iance-covariance matrices of the observation sets. 
Infinitely large variances are applied to non-weighted 
parameters resulting in zeros in the corresponding diagonal 
elements of P2. 

After minimizing 4 > and eliminating the unknowns Ki, K2, 
Vj, and V2, the least squares estimate for X is given by the 
solution of the normal equations 



(Ai^PjAi + A2‘fp2A2)X + (AjPiLi + A2^P2L2> 



0 



Since, however, the A2 matrix arises from direct observa- 
tions on the parameters, the partial derivatives with 

OF2 

respect to the parameters = I, the identity 

matrix and thus the above equation reduces to the form 



(Ai^PiAi + P2)X + (AjPiLi + P 2 L 2 ) = 0 

This has the solution 

A 



X = -(A/PiAi + P 2 )‘‘ + (Ai^PjLi + P2L2) 



The variance covariance matrix of the parameter estimates is 
obtained by normal error propagation methods and is given by 

\ = (A/PiAi + P2)-i 

with the a posteriori variance of unit weight by 



A 

Oo 



2 






"1 



"2 



in which ni equals the number of [signal time interval] 
observations, n2 equals the number of a priori parameter 
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observations, and u equals the number of parameters in the 
adjustment. . . . 

If the assumptions going into the adjustment are good, 
then the a posteriori variance should be close to 1.0. This 

statistic is an indicator of the quality of the adjustment. 

The least squares adjustment technique is an excellent 
tool for determining solutions of unknown variables. However, 
certain limitations do apply. The process assumes that all 
systematic errors have been removed or accounted for, and that 
all remaining errors are randomly distributed. Systematic 
errors that still unaccountably exist in the observations will 
bias the adjustment such that it may attempt to distribute the 
errors to all the observations and shift the parameters. 

Other factors which may degrade the adjustment are: 

* A poor physical model such as one that ignores scale 
error. 

* The incorrect weighting of observations. 

* Small residuals which may be a result of poor network 
geometry or insufficient observations rather than good 
observations . 

For additional information on adjustment computations, see 
Uotila (1986) for derivations of appropriate expressions for 
least squares adjustments which use a variety of different 
systems or groups of observations with their associated 
constraints . 

Appendix A gives greater detail on the least squares 
adjustment process for this Pegasus data analysis. 



35 



V. COMPUTATIONAL PROCEDURE 



The data set used for most of the results below was the 
downcast data from Pegasus Drop #157 on August 3, 1989. 
Occasional comparison studies were also made with the #157 
upcast data, and that of #158 which occurred in the same 
vicinity later the same day. 

As stated in Chapter II, Pegasus data is recorded every 16 
seconds as the instrument descends and ascends during a 
deployment. For this study, the portion of the depth profile 
of Drop #157 from about 16 m to 2200 m was used, providing 306 
Pegasus records. This carried the instrument from surface 
waters down through the plane of the transponders. 

It needs to be stressed at the outset that the 
computational procedures adopted in this thesis should not be 
considered as optimum, but rather were developed as processing 
proceeded in order to overcome difficulties encountered with 
the specific data set associated with the CIO array. If the 
problems ultimately discovered with the data set had been 
known at the beginning, then some procedural points in the 
data processing would have been slightly revised. This aside, 
the purpose of the study was to demonstrate the capabilities 
of least squares estimation procedures to resolve both unknown 
transponder positions and Pegasus velocity components. As 
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will be seen in the following discussions, this was more than 
adequately achieved. 

Because of the poor initial positions of the transponders, 
it was decided to select a small data set of 74 Pegasus 
records from the total of 306 collected during drop #157 and 
use these to help provide improved transponder coordinates. 
This data set was determined by dividing the full depth of the 
Pegasus drop into 36 intervals and taking a pair of 
consecutive records in each interval. Using the a priori 

standard deviations derived in Chapter III, an initial 

solution for the 74 Pegasus positions and the four 

transponders was obtained. The very small positional standard 

deviation associated with T1 essentially served to hold this 

transponder fixed in the resulting solution. 

As described in Chapter IV, (Ai^P^Ai + P 2 ) is the variance- 

covariance matrix of the adjusted parameters. The upper 

tridiagonal portion of this matrix is stored in vector form 

and can be accessed to give information on the variances and 

covariances of the adjusted parameters as shown in Figure 13. 

For example, the portion of this figure which describes 

Pegasus Position 1 gives the variances of the x, y, and 2 

coordinates, , and Op^j^, respectively. The other 

three elements in this upper tridiagonal, Op p , Op p , and 

xl yl xl zl 

Op p , provide the covariances between the x, y, and z 
yl zl 

coordinates of the first Pegasus position. Similarly, the 
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Figure 13. Variance-Covariance Matrix 
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covariances between the coordinates of different positions are 
provided in their appropriate locations as shown. The 
variances and covariances of the transponder coordinates are 
located in the lower right hand corner, and among themselves 
comprise an upper tridiagonal 12 x 12 submatrix. 

The adjustment solution is also used to compute observa- 
tional residuals which in turn are used to compute an a posteriori 

variance of unit weight. For the 74 record run this a posteriori 

variance turned out to be 0.1196, a very small number compared 
to that which should ideally have been close to unity. This 
a posteriori variance when multiplied through the variance 

covariance matrix enables the matrix to be scaled in order 
that it provide supposedly unbiased estimates of the parameter 
variances and covariances. This was done with this data set 
and the resulting transponder coordinates with their newly 
estimated standard deviations (of the order of 10-20 m) taken 
and used as a priori information in the final solution in which 

all 306 records were used simultaneously. It was felt that 
this procedure would enable the transponder coordinates with 
their associated accuracy estimates to more closely represent 
the type of situation usually found in a good transponder 
network. 

In retrospect, however, it appears that it may have been 
more appropriate to have run the 74 point data set with = 
0.0005 s rather than the 0.002 s actually used. When, toward 
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the very end of this study, a Ofime of 0.0005 s was allocated to 
the measured time delays in the 74 point data set, the a posteriori 

variance of unit weight became 0.96, and the resulting 
standard deviations on the transponder coordinates 
approximately 40 m. The coordinate solutions for the 
transponders did not change by more than 5 m in position 
although the computed depth of T4 did increase by a further 
10 m. 

Although this second solution is to be preferred over the 
first, it was felt that the work and time involved in altering 
all the results and documentation already completed was not 
justified, given the minimal impact that it would have on the 
final results. In fact, it was clear that it would not change 
any of the conclusions resulting from this study. 

Ideally, when dealing with data from a number of Pegasus 
drops on a single transponder network, it would be best to use 
a small, but representative (in depth) data set from each 
drop, and merge these together in a single adjustment to 
provide an optimum set of transponder coordinates for that 
array which could, where appropriate, be held fixed in 
subsequent simultaneous processing of complete drops. 
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VI . RESULTS OF STATION CIO ADJUSTMENT 



A. TRANSPONDER COORDINATES 

With improved transponder coordinates and their associated 
standard deviations determined through the 74 record procedure 
described in the last chapter, the entire 306 record data set 
was then run. The results from this full analysis are 
discussed below. 

The final transponder positions were moved (in total) from 
their original Loran estimates as shown in Figure 14 and Table 
2 below. 
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Figure 14 . Transponder Movement 
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TABLE 2 



MOVEMENT OF TRANSPONDER POSITIONS 



Transponder 


X 


Y 


Z 


T1 


0 


0 


4.5 


T2 


-20.3 


-41.9 


-6.2 


T3 


32.2 


116.0 


-7.1 


T4 


-90.8 


149.3 


-6.5 


; must be stressed, however, that 


this data 


proved to be 


internally consistent. 


Using 


the same 


transponder 



positions (as determined by the 74 Pegasus position data set) 
to adjust data from another drop #158 (which used only 13 
Pegasus positions) , physically nearby and later the same day, 
produced somewhat different transponder coordinates although 
it converged within itself. The resulting Table 3 is shown 
below. 

It is suspected that the reason for this inconsistency 
lies both in the very low degrees of freedom in the adjustment 
of drop #158 (leading to a statistically weak solution) and in 
the overall weak a priori positions of the transponders. In a 

well-surveyed four transponder array, this problem would not 
exist. 

It appears certain that the location of Pegasus within 
this weak transponder array has an effect on the final 
adjusted transponder positions. With the origin held fixed in 
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TABLE 3 





TRANSPONDER 


COORDINATE 


SOLUTIONS 






Trans- 

ponder 


X 


Y 


Z 




T1 


0 


0 


1882.86 




T2 


-1387.53 


1530.88 


1987.15 


Drop a 157 


T3 


-1769 . 38 


-66.52 


2229.11 




T4 


-540.12 


2049 .71 


1893.48 




T1 


0 


0 


1879.13 




T2 


-1244.88 


1502.60 


1991.22 


Drop #158 


T3 


-1779.68 


-181.94 


2232 . 52 




T4 


-533 . 65 


1950.66 


1896.60 




T1 


0 


0 


-3.73 


Difference 


T2 


-142 . 65 


28.28 


-4 . 07 


Between 

#157-#158 


T3 


10.3 


115.42 


-3.41 




T4 


-6.47 


99.05 


-3 . 12 



both cases, the array tended to skew in the direction of the 
drop (see Figure 15) . A set of coordinates which would be 
appropriate for all drops might be obtained by making a series 
of Pegasus drops out along the edges of the array near the 
centers of the baselines, in addition to the drop in the 
center. 
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Figure 15. Adjusted Array from #157 and #158 Drops 

It is well to note, however, how the transponder positions 
have been improved over those of our originally assumed 
coordinates, especially those of T4 whose position was largely 
unknown (see Table 4) . 

The horizontal drms values were (when using = 0.0005s): 

* T1 = 0 (held fixed) . 

* T2 = 11.2 m. 

* T3 = 13.2 m. 

* T4 = 12.4 m. 

These drms values suggest an improved horizontal precision 
of 13.5 meters or less. 
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TABLE 4 



STANDARD DEVIATION OF TRANSPONDER POSITIONS 





Transponder 


X 


Y 


Z 




T1 


0 


0 


10 


Initial 


T2 


100 


100 


10 


Standard 

Deviations 


T3 


100 


100 


10 




T4 


100 


100 


100 




T1 


0.01 


0.01 


1.88 


Standard 


T2 


8 .32 


7 .48 


1.75 


Deviations 

After 


T3 


6.79 


11.27 


2 . 03 


Adjustment 


T4 


11.67 


4 . 12 


2.78 



B. PEGASUS POSITIONS 

The precision of the Pegasus positions is determined from 
the variance-covariance matrix as described in Chapters IV and 
V. Typical values at various depths are shown in Table 5. 

If geometric studies are desired, these variance- 
covariance values produce a tri-axial error ellipse for each 
Pegasus position. For ease of perception, the axes of these 
ellipses can be converted to an orthogonal system through the 
use of eigenvalues and eigenvectors. (For a discussion of 
this process, see Mikhail, 1976.) 
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TABLE 5 



PRECISION OF PEGASUS POSITIONS 



Depth (m) 


Ox 


Oy 


Oz 


16 


4.92 


6.32 


0.51 


248 


4.61 


6.18 


0.49 


503 


4.31 


6.00 


0.50 


750 


3.96 


5.80 


0.53 


1000 


3.70 


5.62 


0.58 


1250 


3.44 


5.56 


0.67 


1502 


3.31 


5.52 


0.86 


1750 


3.24 


5.16 


1.29 


2005 


3.53 


5.06 


1.98 


2183 


3.80 


4.76 


1.53 



C. CURRENT VELOCITIES 

Pegasus velocities were determined by using the adjusted 
values of the Pegasus positions which were determined by the 
methods described in Chapter V. 

It is well to note here the similarities and differences 
between the Pegasus positions as determined by this 
adjustment, and those determined by the method currently in 
use by the NPS Oceanography Department. 

The initial Pegasus positions used to start this adjust- 
ment had been derived by combining three observation equations 
into two, and then solving the two equations for the unknown 
X and y of the Pegasus position. (See procedure explained in 
Appendix A.) These positions were then refined by the 
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adjustment process using four transponders where Pegasus depth 
was constrained by pressure. 

The present method used by NPS as described in Chapter II, 
uses only two observation equations to solve for the positions 
(i.e., using only two transponders), without further 
adjustment and without the pressure/depth constraint. 
Velocities derived from these two sets of positions are 
compared in Figures 16 and 17. (Both sets of velocities have 
been filtered by a seven point running average.) Figure 16 
plots the U components of adjustment derived velocities 
against depths, those derived by using the present NPS method, 
and finally the differences between the two approaches. 
Similarly, Figure 17 plots comparable information for the V 
components . 

The profiles of the two methods exhibit very similar 
configurations through much of the range, with notable 
exceptions occurring at the surface, at mid-range (see 
discussion of residuals in the section on depth/pressure 
analysis which follows) , and at the bottom of the cast. 

The differences in surface velocities between the two 
methods may be a result of using different sound velocity 
profiles. The large divergences at depth where the adjusted 
speeds peak at about 22 cm/s and those of NPS at 111 cm/s, are 
almost certainly due to a loss of precision in the two- 
transponder solution, as well as inconsistencies in the 
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Figure 16. U Component 
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Figure 17 . V Component 
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original pressure measurements. Table 6 provides a sample of 
the velocities derived by the two approaches, and shows the 
differences between them. These differences begin to increase 
markedly at about 2100 m, indicating the substantial role the 
adjusted transponder network and the inclusion of a pressure 
observation play in the solution. 



TABLE 6 

VELOCITY COMPARISONS 



Depth 
(Approx. ) 
(m) 


Adjustment 

derived 

velocity* 

(cm/s) 

(1) 


NPS 

Velocity 
( cm/s) 
(2) 


Velocity 
Dif fer- 
rence 
(cm/s) 
(2)-(l) 


Differ- 
ence 
in U 
(cm/s) 


Differ- 
ence 
in V 
(cm/s) 


40 


14.23 


13.24 


-0.99 


+ 6.30 


+ 5.49 


112 


18.89 


25.41 


+ 6.52 


+7 . 09 


+ 6.29 


660 


2 . 10 


1.65 


-0.45 


-0.35 


+ 0.29 


1280 


10.01 


10.74 


+ 0.74 


-1.01 


-0.77 


1380 


1.08 


9.45 


+8.37 


-7.06 


-7.74 


1820 


10.34 


12.66 


+2 .33 


+ 1.93 


+ 1.30 


2090 


17.32 


19.49 


+2.16 


+ 0.81 


+2.20 


2112 


19.58 


35.84 


+16.27 


+11.61 


+11.44 


2145 


23 . 07 


68.05 


+44.98 


+33.89 


+30. 10 


2175 


16.36 


110.75 


+94 . 39 


+71.56 


+62.95 


2190 


16.36 


49.76 


+33.40 


+25.57 


+22 . 62 


♦Adjustment using o 


t = 0.002, 


o = 2 . 






Seven 


point filter 


used for 


both sets 


of velocities. 
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D. PEGASUS VELOCITY ERRORS 



Estimated velocity errors were obtained through the 
propagation of positional errors found in the variance- 
covariance matrix of the final run, yielding o^, o^, and o„, the 
errors in the x, y, and z directions, respectively. This 
propagation of errors proceeds as follows. 

The desired velocity components are given by the 
expression: 





u 




(X2-xJ/At 


= 


V 


= 


(Y2-yi)/At 




w 




(Z2-Zi)/At 



where (Xi,yi,Zi) and are the positions of Pegasus at 

respective 16 second intervals (i.e.. At = 16 s) . 

The variance-covariance of V is given by 



— 










Ouv 


^UW 




^VU 




^ vw 


— GS^pG^ 


^wu 








— 









where 



d/dx^ 3/3yj_ d/dz^ 



3 / 8 X 2 3/3y2 3/3Z2 



G = 




0 

-1 

0 



0 

0 

-1 



1 

0 

0 



0 

1 

0 



0 

0 

1 
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and 




1P2 

In this last expression EP^^ and EP^ are the variance- 
covariance matrices for Pegasus positions 1 and 2, 
respectively, while ^ 1^2 ^ ^2^1 their cross 
covariances. These are obtained from the full adjustment 
variance-covariance matrix described in Chapter V. (Appendix 
B contains an algorithm for computing the variance-covariance 
matrix of the current velocities from data obtained from the 
subroutine MAT in the Least Squares Adjustment program.) 

These values provide information on the precision of 
Pegasus velocity at each position, and change according to the 
geometry between Pegasus and the transponders. Table 7 
provides velocities and their standard deviations at various 
depths as computed by the adjustment. The horizontal velocity 
errors are greatest near the surface, and improve as Pegasus 
descends toward the transponders where, in the horizontal 
sense, the geometry becomes tighter. For the vertical 
velocity, the reverse is true. 

These velocity errors were disappointing. At the surface, 
a velocity in the 'X direction of -17.8 cm/s had a standard 
error of ±10.7 cm/s. In the Y direction, a velocity of 



52 



TABLE 7 



VELOCITIES AND STANDARD DEVIATIONS 



Av. 



Depth 

(m) 


U (m/s) 


V (m/s) 


W (m/s) 


Ou (m/s) 


Ov (m/s) 


o„ (m/s) 


20 


-0.1784 


0.0169 


0.4347 


0.1073 


0.0822 


0.0403 


230 


-0.0886 


0.0789 


0.4566 


0.0988 


0.0759 


0.0408 


740 


0.0236 


0.0204 


0.4511 


0.0795 


0.0621 


0.045 


1233 


-0.0729 


0.0011 


0.4141 


0.0643 


0.0511 


0.0573 


1739 


0.1178 


-0.0238 


0.4734 


0.0578 


0.0439 


0.1107 


2180 


0.1989 


0.0681 


0.3379 


0.0569 


0.0454 


0.1339 


1.7 cm/s had a 


standard 


error of 


±8.2 cm/s. At 


the bottom 



results were somewhat improved with u = 19.9 cm/s and = 5.7 
cm/s, V = 6.8 cm/s and = 4.5 cm/s, but not substantially so. 

These velocities and error estimates were obtained by 
looking at only individual pairs of positions. The NFS 
Oceanography Department computes velocities by using a seven- 
point running average for each Pegasus recorded depth. This 
same filtering technique used with the adjusted data would 
improve precision by l/'H^ or 0.38, if positions were 
uncorrelated. However, these positions are correlated, so the 
expected improvement is somewhat smaller. 

There are two ways the precision of the unfiltered 
velocities could be improved: 

* Use longer time intervals, so that as At becomes larger 
in the error propagation equation discussed above, the 
variance-covariance matrix 2 V diminishes. 
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* Average the 16-second velocities (in the same manner as 
NPS) by using a seven-point running average which would 
improve the precision if the velocities were 
uncorrelated . 

However, both of these techniques, while gaining precision, 
sacrifice vertical resolution. 

There are two important aspects to these formal velocity 
errors. In the first instance the signal time intervals are 
only measured to 0.0001 s. At a conventional sound velocity 
of 1500 m/s this is equivalent to 15 cm in the derived range. 
This in turn propagates into a velocity error of approximately 
1 cm/sec between consecutive Pegasus positions if the 
transponder positions are assumed to be without error. 

In the second instance, the formal errors in the 
transponder positions propagate directly into errors in the 
Pegasus positions and thence into the velocity components. 
These latter errors have by far the most significant influence 
on the derived standard deviations of the velocity components 
for Pegasus. This in turn provides an added emphasis on the 
need for a strong transponder survey. 

E. DEPTH/ PRESSURE ANALYSIS 

To determine at what depth the pressure/depth measurement 
becomes critical in solving for current velocities, a separate 
computer run was made with the standard deviation of the 
Pegasus Z coordinate changed from 2.00 to 50.00 m. The 
results are graphically illustrated in Figure 18. 
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SPEED (CM/S) 




Here, the U and V velocity components derived from the 

adjustment using Op = 2 m are plotted on the left; speed 

z 

differences between the U and V components of the two sets of 

adjustments using Op =2, and Op = 50, respectively, are 

z z 

illustrated on the right. The differences are minimal down 
through the water column until about 1700 m as Pegasus 
approaches the plane of the transponder. DU reaches a maximum 
positive difference of 1.16 cm/s at 2008 m, and subsequent 
maximum negative difference of -3.53 cm/s at 2135 m. 
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Respective values for DV are 0.49 cm/s at 2001 m and 
-2.26 cm/s at 2135 m. 

An analysis comparing the standard deviations of Pegasus 
depth positions shows a large maximum of 36 m at a depth of 
2000 m as derived from the free-floating adjustment, 
compared with a maximum of 2 m at a depth of around 1940 m for 
tightly-held P^. The latter reflects the fact that the a priori 

precision estimate for pressure observations was itself 2 m, 
and that travel time measurement geometry gives no additional 
information on the adjusted z value. 

Table 8 shows observational residuals after both of the 
adjustments and compares the two sets of residuals at similar 
depths. (These figures compare runs made with Oii„e = 0.002 s, 
rather than the final 0.0005 s.) 

At the surface both sets of residuals are small. Where 
the pressure/depth relationship is tightly constrained, 
residuals were highest at the bottom of the cast where Pegasus 
moved through the plane of the transponders, i.e., between 
about 1870 m to 2240 m. Here the resolution of depth becomes 
less certain due to the poorer geometry of the Pegasus- 
Transponder array. These large residuals at depth indicate 
there may be a bias in the pressure measurements, thus 
suggesting the presence of a systematic error in the pressure 
sensor (perhaps due to a hysteresis in temperature readings) 
that is unaccounted for in the adjustment model. 
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TABLE 8 



OBSERVATIONAL RESIDUALS 



z 




T1 


T2 


T3 


T4 


Surface 


2 


-0.0001 


-0 . 0002 


• 0.0001 


0.0001 






0.0000 


0 . 0000 


0.0000 


0.0000 






0.0000 


0.0001 


-0.0001 


-0.0001 




50 


0.0000 


-0.0002 


0.0001 


0.0001 






0.0000 


0.0000 


0.0000 


0.0000 






0.0000 


0.0001 


0.0000 


0.0000 


1400 m 


2 


-0.0004 


-0.0012 


0.0006 


0.0011 


(approx. ) 




-0.0005 


-0.0013 


0.0006 


0.0012 






-0.0005 


-0.0012 


0.0006 


0.0011 




50 


-0.0004 


-0.0011 


0.0006 


0.0010 






-0.0004 


-0.0012 


0.0006 


0.0012 






-0.0004 


-0.0013 


0.0007 


0.0012 


Bottom 


2 


0.0014 


0.0011 


0.0007 


0.0003 






0.0016 


0.0014 


0.0008 


0.0004 






0.0019 


0.0016 


0.0008 


0.0003 




50 


0.0001 


0.0001 


0 . 0000 


0.0000 






0.0002 


0.0002 


0.0000 


-0.0001 






0.0002 


0 . 0003 


-0 .0001 


-0.0002 



Residuals at the bottom of the cast from the adjustment 

where depth was allowed a very loose constraint are uniformly 

low. The differences between these two sets of residuals at 

the bottom, as well as differences in velocities and an 

increase in Op at this depth as discussed above, all tend to 
z 

corroborate the suggestion of a bias in the pressure 
measurements . 
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Residuals also increased in both cases for a few positions 
around the depth of 1400 meters. This appeared to occur where 
the trajectory of the instrument reached its most westerly 
position. Reasons for this are unclear, although a possible 
explanation may be related to the fact that this depth and 
position approximate the location of the shelf edge of the 
canyon. 

It is also possible that one or more of the travel time 
observations in this part of the drop had large errors. This 
interpretation is further suggested by the anomalous current 
shear computed at this depth from the original data (see 
Figures 16 and 17) . 

It should additionally be noted that these large residuals 
at depth and at mid-range seem to display a systematic 
structure wherein their sign is consistently the same. If the 
errors were truly random, the positive and negative signs of 
the residuals would be distributed more or less evenly. This 
situation is another indication of some systematic error which 
has not been accounted for in the adjustment model. 
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VII. CONCLUSIONS AND RECOMMENDATIONS 



It is unfortunate that the data set and the transponder 
array used in this study lacked the quality desirable for de- 
tailed analysis. Many of these difficulties were, unfortu- 
nately, only discovered as data processing proceeded. It is 
encouraging to note, however, that the least squares proce- 
dures described here enabled the identification of problems 
with the data which would otherwise have passed undetected. 

From a theoretical standpoint, the least squares process 
must provide better results than the standard NFS positioning 
techniques because of the fact that it uses all the observa- 
tional data. It has the added advantage of providing vital 
statistical information on the quality of the derived results. 
Unfortunately, in the case of the CIO network, the uncertainty 
regarding the original positions of the transponders made it 
difficult to perform a comprehensive comparison between the 
least squares technique and the standard NFS methodology. 

It is observed, however, that the least squares method 
obtained a solution for transponder positions that converged 
to less than one meter, with standard deviations that were 
much smaller than those with which the adjustment started. 
Precision of the transponder coordinates showed horizontal 
drms values of less than 15 meters, with standard deviations 
of transponder depth less than three meters. While these 



59 



figures may be optimistic (see Chapter IV) , they are 
indicative of the improvements which can be achieved using 
these methods. 

Regarding the precision of the horizontal Pegasus 
positions, the standard deviation of the X values ranged from 
4.9 down to 3.2 m, reaching a minimum around a depth of 1800 m 
just at the upper limit of the transponder plane, and then 
increasing slightly to the bottom of the cast. Pegasus Y 
values showed standard deviations of 6.3 to 4.8 m, with 
passage through the transponder plane not as noticeable. 

Standard deviations of the Pegasus Z values, when they 
were tied tightly to pressure, ranged from 1.4 m near the 
surface to a maximum of 2.0 m at a depth of 1943 m. When the 
depth was only loosely constrained, standard deviations varied 
from 2.3 m near the surface, to 8.6 around depths of 1800 m 
just prior to entering the plane of the transponders, and 
reached a maximum of 36.5 at a depth of 2000 m. From there 
they steadily diminished back to 10.4 at the final depth of 
2165 m . 

Comparison of the results of the NPS method and the least 
squares method shows a considerable difference in velocities 
at depths where Pegasus passes through the plane of the 
transponders. This difference at depth occurs also in 
velocity comparisons between adjustments made holding depth 
constrained with a standard deviation of 2.0 m, and allowing 
it to float more freely with a standard deviation of 50.0 m. 
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This leads to a suspicion that there may be a systematic error 
in the pressure reading that has not been accounted for in the 
model, perhaps due to a hysteresis in temperature readings. 
Accuracy of pressure observations becomes critical at depths 
below 1700 m. 

It was disappointing to note that the standard deviations 
on the computed Pegasus velocities (when using At = 16 sec) 
were often of a similar magnitude as the velocities themselves 
(i.e., 5-10 cm/s). However, with resolution of signal travel 
time possible only to the nearest 0.0001 s, resolution of 
range becomes ±15 cm. This in turn propagates into velocity 
errors in the order of 1-2 cm/sec. Therefore, it is unrealis- 
tic to look for much greater precision than that. 

It is recommended that: 

* NPS refine the existing least squares techniques and the 
adjustment software to facilitate its use on a regular 
basis, for production operations. 

* Four transponders be used in each network. This will 
both strengthen the solutions for the Pegasus velocities 
and provide a reasonable measure of redundancy for the 
adjustment procedure. 

* More acceptable positioning of transponders be under- 
taken, within the time constraints involved, including 
observation of all baselines. Close attention should be 
paid to both length and azimuth. 

* That a small, but well-distributed data set of Pegasus 
records be used from each drop to assist in providing a 
unique set of transponder positions for each network and 
that these positions be held fixed in the subsequent 
adjustment of the full set of records from each drop. 

* Where practical, pressure information always be 
collected. It both adds to the overall redundancy in the 
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network and enables additional parameters to be 
introduced if necessary. 

* A very well-controlled set of experimental data be 
collected on a well-positioned five transponder network. 
This will enable clarification of the following issues: 

- The lack of consistency between the adjusted transpon- 
der positions when using different data sets for the 
same area. 

The question of pressure/depth relationships, 
especially at depth. 

- Possible causes for high residuals at a given depth 
(1400 meters in the case of station CIO) . 

The introduction of additional parameters to describe 
possible mis-calibration of the pressure head or 
hysteresis of temperature readings. 

The investigation of the actual motion of Pegasus as 
it descends and ascends. 

The geometric impact of drops made on the edges of the 
transponder array as opposed to those made at the 
center of the array. 
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APPENDIX A 



LEAST SQUARES ADJUSTMENT PROGRAM 

The least squares adjustment equation for a combined 
system of two sets of observations is (Uotila, 1986, p. 97): 



A 

X 



-(Ai^PlAi + A2‘"P2A2)-1 (A/PiLi + A2^P2L2) 



where, in the terminology of ' the program: 



X = column vector of corrections to the parameters 
(NP/1) ; 

L = Observations (NOBS,l); 

A = Jacobian matrix, or the partial derivatives of each 
observation with respect to each parameter 
(HOBS,NP) ; 

a” = Transpose of A (NP,NOBS) ; 

P = V7eight matrix (NOBS, NOBS); 



N = Maximum 


4i 

n 


of 


Pegasus positions to 


be 


determined; 


NT = Maximum 


# 


of 


transponders in the 


network; 


NOBS = Maximum 


i 


of 


observations 






NP = Maximum 


# 


of 


parameters allowed = 


(N 


X 3) + 



(NT X 3) . 
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A. FIRST SET OF OBSERVATIONS 
1 . Li Matrix 

The first set of observations uses one-way travel 
times from the transponders to Pegasus as the observations. 
Parameters are the x, y, and z coordinates of each Pegasus 
position . 



Li = Fi(x) 

Time = Distance/Velocity 

Position 1: 

Timei = [(XPi-Xi)2 + + ( ZP^-ZJ i/y 

Time^ = [(XPi-X2)2 + + (ZPi-Z 2 )"]Vv 

Timea = [(XPi-Xa)^' + (YPi-Ya)^ + (ZPi-Za)2]Vv 

I 

Time^ = [(XPi-XJ^ + (YP^-YJ^ + (ZPi-ZJ^jyv 



Position N: 

TimOj = time intervals from transponders 1,4; 

V = sound velocity; 
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XP,YP,ZP = Pegasus position coordinates; 

Xj,Yi,Zi = transponder coordinates; 

# Observations = N x 4 . 

Each Pegasus position therefore has four time- interval 
observations, one for each transponder. Thirteen Pegasus 
positions will provide 52 time observations and (13 x 3) 
parameters, and therefore 13 degrees of freedom for this first 
set of observation equations. For each additional Pegasus 
position added to the data set, an additional degree of 
freedom is gained. If a priori constraints are introduced for 

the transponder coordinates and the transponder positions 
solved for in the adjustment (as has been done here) , then the 
number of observations and the number of unknowns rises by 12. 

The adjustment procedure requires an a priori knowledge 

of the transponder positions, the depth of each Pegasus 
position (derived from the pressure information) and the sound 
velocity for each Pegasus position. Estimates of the x and y 
coordinates are obtained by taking the three observation 
equations for each Pegasus position that contain information 
from transponders Tl, T2 , T3 (presumably the best known). 
Subtracting the first equation from the second, and the second 
from the third, will yield two equations in two unknowns, and 
therefore a resulting first guess solution for XP and YP at 
that position. 
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The Lj matrix in the adjustment equation is a matrix 
of differences between observed time intervals and what the 
observation equations would yield using our first guess 
parameters . 



Li 



Li 



calculated 



L 



1 observed 



2 . Ai Matrix 

The Ai matrix contains the derivatives of each 
observation with regard to each parameter: 



Ai Matrix 




If 



Ti = [(XPj-Xi)2 + (YPj-Yi)2 + (ZPj-Zi)2]i/v = D/V 

i = 1, 2 . . 4 
j = 1,2 . . N 



Then 



3T, 



axPj^ = 



(XP^-X^) 

DV 
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<ZPn-^4> 

3ZP., DV 

N 

The derivatives of the transponder positions will be 
the negatives of the corresponding Pegasus position 
derivatives . 

3F^ (XPj-X^) 

JiT "" DV 

1 

This will create a matrix in the form: 



Peg . 
Pos 1 



Peg . 
Pos 2 



Peg . 
Pos N 



3/3p^ d/dp^ a/DT^ 





X 


. X , 


• X 














X 


. X 


, X 




















^2 


X , 


. X , 


» X 




















X 


, X 


, X 














^3 


X , 


, X , 


, X 


























X , 


, X , 


, X 








^4 


X , 


, X , 


, X 
































X , 


. X , 


. X 


^ 1 








X , 


» X , 


, X . 








X , 


. X , 


• X 




















^2 








X , 


. X , 


, X . . 














X , 


. X 


. X 














^3 








X , 


, X , 


, X . . 




















X 


. X , 


• X 








^4 








X , 


, X , 


, X . 


























X , 


r X , 


, X 


M 














X , 


r X , 


» X 


X , 


r X , 


, X 




















^2 














X 


rX 


» X 








X , 


, X 


, X 














*^3 














X 


» X , 


, X 














X 


, X , 


, X 






















X , 


. X 


. X 




















X , 


, X , 


, X 
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where : 



ti 


= 


Time^ ; 




Pj 


= 


XPj, YP 


j ' 


Ti 


= 


Xi, Yi, 


Zi 


i 


= 


1,2 . . 


4 ; 


j 


= 


1,2 . . 


N. 



ZP, 



If submatrices are described by 





X, X, X 




x,x,x, 0,0, 0,0, 0,0, 0,0,0 


aj = 


X, X, X 
X, X, X 
X, X, X 


So = 


0,0,0,X,X,X,0,0,0,0,0,0 
0,0,0,0,0,0,x,x,x, 0 , 0,0 
0,0,0,0,0,0,0,0,0,x,x,x 



Then the matrix becomes 



a^ 0 0 ••••••• 0 3^ 

0 a2 0 0 S 2 

0 0 ^3 h ^3 






For each Pegasus position there will be a 4 x 3 block 
for the Pegasus derivatives and a 4 x 12 block for the 
transponder derivatives. All other positions in the matrix 
will be zero. Therefore, to drastically reduce computer 
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memory and time requirements, the program was developed to run 
in a sparse matrix form whereby only the non-zero elements are 
used. To do this, the whole A matrix is divided into parts: 

* A — which contains the values for the Pegasus components. 

* S--which contains those for the transponders. 

3 . P , Matrix 

Pj, the weight matrix, is set up as a diagonal matrix 
and stored in column vector form, assuming no covariances 
between observations. If it is desirable to add covariances, 
then an error propagation subroutine can be inserted to fill 
out the weight matrix and the program altered accordingly. In 
our case the chosen standard deviation for the time interval 
was 0.002 s, later changed to 0.0005 s. It was entered into 
the program in the error propagation section where the code 
reads: 

Do 35 I = 1,NT 

P(l,l) = 1.0D0/0.0005DO**2 
35 Continue 

B. SECOND SET OF OBSERVATIONS 
1 . L, Matrix 

The second set of observations corresponds to various 
parameters which have been observed, in this case the position 
coordinates of the four transponders, and the position of the 
Z coordinate of each Pegasus position which allows the depth/ 
pressure relationship to be constrained. 
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^2 — F" ( ^ ) 2 

Xi = Xi 

Yi = Yi 

. X,Y,Z of each transponder 

Z, = Z, 

ZPi = ZP^ Depth of Pegasus as computed 

i = 1,N from pressure. 



As described above, we already have first guess 
estimates of all these observations. The L 2 matrix will 
initially be zero since for the first solution , the observa- 
tions equal the parameters, i.e., Xj = Xj. However, with 
further iterations, this will not be the case. 

2 . A, Matrix 

The A 2 matrix contains the derivatives of each 
observation with regard to each parameter: 



A 2 



= ^ 



3X 



a 



For example: 




1 



8x 

"^11 



1 

else 



0 
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This produces the following diagonal matrix of I's and 



0 ' s: 



0 

0 



1 

0 

0 



1 




1 

1 



3 . P-, Matrix 

P 2 is a similar diagonal matrix stored in column vector 
form. It weights all Pegasus z positions to allow for 
tightening or loosening the depth/pressure relationship. This 
chosen standard deviation, in our case a value of 2.00, enters 
the program through an interactive request at the beginning of 
the program run. The standard deviation of each transponder 
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position coordinate is road into the program from the original 
data set. The matrix will appear as: 



1 / 2 

zl 

U 



V. 2 



z 2 



•4 2 



'xl 



1/ 2 
' o 1 
yi 



V. 2 






C. MATRIX OPERATIONS 

Now all the matrices are in place and matrix operations 
can begin. 

X = -(Ai’'PiAi + A2’'P2A2) (A/PjLi + A2’'P2L2) 
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since the A 2 matrix is a diagonal of O's and I's, the 
second term in the first parenthesis merely adds weight to the 
appropriate diagonal position in the first term. 



D. REQUIREMENTS FOR RUNNING THE PROGRAM 
1 . Initial Steps 

Before the program is run, the following steps must be 



taken: 

* Variables must be dimensioned adequately, following the 
definitions clearly stated in the preface to the program. 

* The value of must be defined (in our case 0.0005), 

and/or the appropriate error propagation subroutine 
added . 

* A data set file must be available for access and set up 
in the following fashion; 

- Transponder coordinates, each with its own standard 
deviation. 

For each Pegasus position: pressure, depth, four 

transponder one-way time intervals, and sound 
velocity. (Sample data set is shown below.) 
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Transponder 

Coordinates 


T1 

T2 

T3 

T/i 


0 

-1324 

-IbOl 

-449 


.0 
.6 
. 6 
. 3 


0.0 

1485.8 
-182. S 
1900.4 


1078. 
1 993 . 
2236. 
1900. 


4 

4 

o 

0 
















Transponder 


Tl 


00 


.01 00.01 




10.00 
















Standard 


T2 


50 


.00 50.00 


10.00 
















Deviations 


T3 


50 


.00 50.00 


10.00 


















T4 


50 


.00 50. OC 


1 


100.00 
















Pegasus 


1 


30. 


3 


30.6 


1 . 


4145 


1 . 


5486 


1 . 


6277 


1 . 


6803 


1485. 


310 


Positions 


2 


2i3 . 


4 


211.7 


1 . 


3255 


1 . 


4245 


1 . 


5191 


1 . 


,5704 


1485. 


100 




3 


593 . 


0 


394.7 


1 . 


2356 


1 . 


3083 


1 . 


4105 


1 . 


4712 


1485. 


140 






SOI . 


5 


576.4 


1 . 


, 1414 


1 . 


.2047 


1 . 


, 3062 


1 . 


,3844 


1405. 


,250 




, 


764 . 


7 


757.7 


1 . 


.0457 


1 . 


.1140 


1 . 


.2043 


1 , 


.3106 


1485, 


,650 




, 


944 . 


0 


934.9 


0. 


.9592 


1 . 


.0202 


1 . 


.1121 


1 , 


.2403 


1486 . 


.200 






1 12S. 


0 


1114.5 


0. 


.8941 


0. 


.9381 


1 . 


.0179 


1 . 


. 1 754 


1486, 


.820 






1 • C6 . 


1 


1292.4 


0, 


.8405 


0 


.8b3l 


0 


.9243 


1 


. 1302 


1487, 


.570 






ues . 


7 


1472.5 


0 


. 7990 


0 


.7 944 


0 


.8472 


1 


.0851 


1489 


.440 






1 667. 


0 


1648.1 


0 


.7620 


0 


. 7455 


0 


. 7922 


1 


.0491 


1489 


. 370 






1846 . 




1824.5 


0 


.734 1 


0 


.7229 


0 


.7534 


1 


.0306 


1490 


.400 






202 1 . 


9 


1997.4 


0 


. 7337 


0 


.69 73 


0 


.7584 


0 


.9940 


1491 


.510 






2200 . 


7 


2173.1 


0 


.7415 


0 


. 6950 


0 


.7910 


0 


.9671 


1492 


.760 






Press- 

ure 


Depth 




^2 




^4 


Sound 

Velocity 



Time Intervals 



* FILEDEFS for input and output files must be in place, 
either typed in before the program starts or available 
through a program exec. The program will run in either 
WF77 or F0RTVS2 fortran. 

- FILEDEF 01 DISK fn ft fm (i.e., Test Data A). 

- FILEDEF 02 DISK(recfm vb Irecl 132 blksize 134). 

- (For large data sets, use FILEDEF 02 DISK(recfm fb 
Irecl 132 blksize 13200) . 

* Be sure enough memory is available. Two M sufficed for 
running at least up through 74 Pegasus positions, but 
4096 K were required for the 306 position run. 

The program starts off by asking for operator input 

(sample answers in parentheses) , requesting: 
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* Project or run description. 

* Number of transponders (4) . 

* Global standard deviation to be given to the Pegasus z 
coordinates (2.00). 

* Convergence limit for adjustment and number of iterations 

(entered separately) (1.000), (4). 

From then on the program works on its own, finally 
producing whatever output has been requested in the output 
file. 

2 . Brief Program Outline 

* Matrices are zeroed out at appropriate times. 

* Pegasus data set is read in. 

* Subroutine APPRO is called to calculate initial approxi- 
mate coordinates of Pegasus. 

* The diagonal weight matrix P^ is set up and stored in 
matrix BIGP. 

* Program constants are calculated. 

* Weight matrix P 2 is generated. 

* The Lj matrix is formed. 

* The matrices A and S are formed by calling subroutine 
SETUP. 

* The normal matrix equations are formed block by block by 
calling subroutine NORM. Outputs include the normal 
matrix (ATPA) , in column vector form called anorm, and 
xhat the ATPL vector called xhat (which subsequently 
becomes the solution vector) . 

* Subroutine P2L2 is called to compute the P 2 L 2 matrix and 
add it to the column vector. It also increments the 
diagonal elements of the normal matrix to include the 
influence of P2. 

* The normal equations are solved by calling a user 
supplied routine that will accomplish the required matrix 
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operations for producing the adjusted values with which 
to correct the parameters. The IMSL routine LINV3P was 
used in this study. Utilizing a symmetric storage mode, 
this algorithm replaces the ATPA matrix by its inverse 
which is the variance-covariance matrix of the Pegasus 
and transponder positions. This (ATPA)'^ matrix can be 
written to file for later access in order to produce the 
error matrix required for current velocity analysis (see 
Chapter V) as used in the program provided in Appendix B. 

* The parameter corrections are tested to see if iteration 
is required. If so, the program then updates the 
parameters and goes through the process again for as many 
iterations as are called for, or until the convergence 
limit has been reached. 

* Residuals and the a posteriori variance of unit weight are 
computed by calling subroutine RESID. (Residuals here 
are, of course, given in units of seconds since they 
relate to time interval observations.) 

3 . Program Outputs 
Outputs include: 

* Original data set and interactive information given at 
the beginning of the program. 

* Adjusted parameters following each interaction. 

* Residuals and a posteriori variance of unit weight. 

* Any other information requested to be written to a file. 

For this study a piece of code. Subroutine MAT, was 
attached at the end of the program to pick out from the anorm 
vector only the data required for current velocity analysis. 
This subroutine is called after the write statement for the a 

posteriori variance of unit weight which uses format statement 

#175. The anorm vector is the upper tridiagonal of the 
(ATPA) variance-covariance matrix stored in column vector 
form. (See the diagram in Chapter V.) 
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The position of the variance of each parameter can be 



obtained by: 



N X (N +1) 

2 

where N is the place number of the parameter. For example, 
suppose the variance for the Pegasus y coordinate in the 
second record was requested. YP 2 is the 5th parameter — XPj, 
YPi, ZPi, XP 2 , YP 2 , (5*6)/2 = 15. The variance for YP 2 , 

therefore, will occupy the 15th place in the anorm vector. 

XPj YPi ZPi XP2 YP2 ZP2 

1 2 

Peg. Pos. 1 3 

Peg . Pos . 2 

The subroutine Mat picks up only those values in the 
anorm vector that fill the 3x3 variance-covariance matrix of 
each Pegasus parameter, and the 3x3 covariance portion that 
relates each parameter to the next one in line. A sample of 
the output follows: 



4 

5 

6 



7 

8 
9 



10 



11 

12 

13 



16 

17 

18 



19 

20 
21 
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Var-covar 
Pos 1 



Covar 

Pos i-Pos 2 



var-covar 
Pos 2 



Covar 

Pos 2-Pos 3 



Var-covar 
Pos 3 



24 .2<8 7SK5f ilS59<l 18 . 4 5280804 18244140 0 . 9 54 4 I 4 24 9 4 5 1 9 9 02 7 4E- 0 1 

18.4528080418244140 5 9 . 8 92 54 98 5S 7 7 7 7 1 8 4 - 0 . S 9 1 S OS 7 7 7S OSO 1 040S 

0.954 4 1424 9 45 1 9 902 7 4E-01 -0 . S 9 1 S OS 7 7 7S OS 0 1 Q40S 0 . 2 4 0 7 1 4 1 8SS 5288 5 2S4 

22. 7 00 7 92 5 92 788 5S4 4 1 8 . S 9 7 9 1 4S 72 7 7 5 1 24 5 0 . 4 5 1 24 02 5 I 5B4 40S242E-0 1 

18.S4 12829 5 40 10044 9 5 9 . 05 4 4 50 5 1 4 SS 4 78 90. -0.44 4 420407 7 9940 7421 

0.8518874441 9S048052E- 01 -0 . 4 5 9S2 90 4 1 1 8 4 1 888 5 0 0 .S20 7 1 09222054S7 1 29E-01 

24.09997 5 0^bV4 4 5 7842 1 8 . 4 0 96 782 7 1 84 1 0 9 58 0.725 9 1 7 4 4 2 14 784888 SE- 01 

18.4098782718410958 5 9 . 9 1 022 5 7 40S 9 52 5 4 4 -0 . S 9S 09 4 4 4S 1 08 4 5 7 9S 7 

0,72 591 74 421 4 7B4e88SE-01 -0 . S 9S 0 9 4 4 4S I 08 4 5 7 f S 7 0 . 2 4000 90 1 SS 7 5909 5 5 5 

22.S4 71 7S 7 1 95 28 1 941 1 8 . S 70 9 1S8 904S 1 900 1 0 . 48 S 7 7 94 18 74 09S 14 04E-0 1 

18.SS09S07 5S9S484 24 5 9.04949780921 77 121 - 0 . 44S 9S42 1 4 58 5 78 1 2 9 9 

0. 4 54 724 I 0004 18 90 14 7E-0 I -0 , 44 54 9 9422 1 59 40 7S 58 0 . S 1 Q8S746 7 9 5 7 7SS 94 lE-0 1 

25.9448 74 488 488954 7 1 8 . S 94 54 9 904284 4 54 1 0 . S 4 S 9 7 1 25 4 4 7 48 5 78 24E-0 1 

18 .S 94 54 9 90428 4 4 541 S9.90982549S8414441 -0.S9599841895S571994 

0.S4S9 7125 4 4 748 5 7824E-01 -0.S959 96 4189 5$ 5719 94 0 .2S9 1 9 54 907 7 40 75099 



The first set of three lines with three values each is the 
variance-covariance matrix for Pegasus position 1^ which is 
symmetric. The second three lines of three values is the 
covariance matrix for position 1 and position 2 
(nonsymmetric) . The third set of three lines is the variance- 
covariance matrix for Pegasus position 2 (symmetric ) , and so 
on. 

The final 30 lines in the output provide the variances 
of the transponder coordinates and all the covariances between 
them. 
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PROGRAM PEGASUS 



* * 

* COMMENTS GIVING A DESCRIPTION OF THE PROGRAM 

Vr * 

IMPLICIT REAL*8 (A-H.O-Z) 

C 

C @^<a(a(9(3@(^(3[@^@(3(3^(5<3(3(3(^(3^(3(3(3(3^(3@(3@(3(3@0(3@00(3^(3@@(3<3(3(3@^^@@^(9(3@0(3(3@ 



C Q @ 
C @ SET MATRIX VARIABLES TO THEIR MAXIMUM SIZES. THE VARIABLES (§ 
C @ USED ARE AS FOLLOWS; @ 
C @ N = MAX. it OF PEGASUS POSITIONS TO BE DETERMINED (20 ) @ 
C @ NT = MAX. if OF TRANSPONDERS IN THE NETWORK (4) Q 
C @ NP = MAX. if OF PARAMETERS ALLOWED = (N*3) + (N'P’VS) @ 
C @ NOBS= MAX. it OF OBSERVATIONS = N*NT (§ 
C (a NV = MAX. if OF NON-ZERO VALUES IN THE A MATRIX = 6*NT*N @ 
C (3 NN = MAX. it OF ELEMENTS IN THE UPPER TRIANGULAR PORTION (§ 
C (a OF THE NORMAL MATRIX = NP*(NP + l)/2 @ 
C Q @ 
C @ THE FOLLOWING MATRICES MUST BE DIMENSIONED TO THEIR MAX. Q 
C Q SIZE PRIOR TO THE PROGRAM BEING USED. Q 
C @ (§ 
C Q XO(NP), P(NT,NT), BIGP(NT.NOBS) ,LB(NOBS, 1) , TRANS(NT), Q 
C Q X(NT),Y(NT),Z(NT),A(NT,3),S(NT,NT’'^3),VALUE(NV),ATP(NT,3), @ 
C @ ATPA(3,3), ATPS(3,NT’’'3) , ATPL(3,1), STP(NT’'’3 ,NT) ,SV(N) Q 
C (a STPS(NT’'^3,NT’'^3), STPL(NT’'^3 , 1) , L1(N0BS,1), L(NT,1), <§ 
C (a IC0L(NV),IR0W(NV), SPS( NT^’^3 ,NT*3) , SPL(NT^‘3 , 1) , ANORM(NN) , @ 
C (a XHAT(NP),XA(NP),VTP(1,NT),V(NT,1),V1(N0BS),SIG(1,1), <§ 

C @ SDX(NT),SDY(NT),SDZ(NT),P2(NP),L2(NP),LB2(NP),AUX(2) (3 
C @ 0 
C 0 NOTE: THE SAME TRANSPONDERS MUST APPEAR IN EVERY PEGASUS 0 
C 0 POSITION. A RECORD IN WHICH ONE (OR MORE) DROP OUT 0 
C 0 IS NOT ALLOWED. 0 
C 0 0 



c 

DIMENSION XO( 72) , P(4 ,4) , BIGP(4 , 80) ,TRANS(4) ,X(4) , Y(4) , Z(4) , 

1 A(4,3) ,S(4,12) ,VALUE(480 ) ,ATP(4,3) ,ATPA(3,3) ,ATPS(3,12) ,SV(20) , 

2 ATPL(3,1), STP(12,4), STPS(12,12), STPL(12,1), ANORM( 2628), 

3 XHAT( 72),XA( 72) ,VTP( 1 ,4) ,V(4 , 1) , Vl( 80) , SIG( 1 , 1) , SD.X(4), 

4 SDY(4), SDZ(4), P2( 72), AUX(2) 

C 

DOUBLE PRECISION LB( 80,1), Ll( 80,1), L(4,l), L2( 72), LB2( 72) 

C 

CHARACTER*12 PROJ 
C 

INTEGER*2 ICOL(480 ), IROW(480 ) 

C 

COMMON SPS(12,12), SPL(12,1) 

C 

C 0 0 

C 0 READ IN THE FIXED DATA AND ANY ASSOCIATED VARIANCES 0 
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C @ UNITS 5 AND 6 READ AND WRITE FROM AND TO THE SCREEN (3 
C @ FOR INTERACTIVE PROCESSING. UNITS 1 AND 2 READ AND @ 
C @ WRITE FROM AND TO DISK FILES. (3 
C @ (3 



c @<3(3<3<3(3<3(3@a(ci(a<a(a^@(a^@(a@(a@0@@@@(a(a@<a@@@@(a(a(a(a@@@@@(a(a@(a^i^^ 

c 

c 

WRITE(6,*) 'PROJECT OR RUN DESCRIPTION (TO 60 CHARACTERS)' 
READ(l.lO) PROJ 
10 FORMAT(A12) 

WRITE(6,*) 'NUMBER OF TRANSPONDERS (THEIR COORDINATES AND STD. ', 
I 'DEVIATIONS TO BE READ FROM THE DISK FILE)' 

READ(1,*) NT 

READd,*) (X(I).Y(I).Z(I), 1=1, NT) 

READ(1,*) (SDX(I),SDY(I),SDZ(I), I = 1,NT) 

WRITE(6,*) 'THE GLOBAL STD. DEVIATION TO BE GIVEN TO THE ' , 

1 'Z PEGASUS COORDINATES' 

READd,*) SDZP 

WRITE(6,*) 'CONVERGENCE LIMIT FOR ADJUSTMENT AND # OF ITERATIONS' 
READd,*) TOL, NITER 

WRITE(2,15) PR0J,NT,(X(I) ,SDX(I) ,Y(I),SDY(I),Z(I),SDZ(I),I =1,NT) 

15 FORMATC 1' ,///,5X,70('*' ),//,10X,A60, //.5X,70( '*■),//, lOX, 

1 'NUMBER OF TRANSPONDERS =', 12 ,//,' TRANSPONDER COORDINATES AND ', 

2 'THEIR STANDARD DEVIATIONS ',//, 4( lOX, 3(F10. 2 , 2X,F6. 2) ,/),/// ) 
WRITE(2,16) SDZP 

16 FORMATC '.SX.'THE GLOBAL STANDARD DEVIATION FOR THE PEGASUS ', 

1 'Z COORDINATES =' ,F7. 2) 

WRITE(2,17) TOL, NITER 

17 FORMATC ',//, SX.'THE CONVERGENCE LIMIT FOR THE ADJUSTMENT =' , 

1 F7. 3,//,6X, 'THE NUMBER OF ITERATIONS ALLOWED =' , 13,///, 

2 SOX, 'INPUT DATA' ,//,5X,' PRESS. ' ,4X,'ZP TIME DELAY l',2X, 

3 'TIME DELAY 2',2X,'TIME DELAY 3',2X,'TIME DELAY 4', 5X.' SOUND ', 

4 'VELOCITY' ,2X, 'APPROX. PEGASUS POSITIONS(X,Y,Z) ' ,/) 

C 

C @ @ 

C @ READ IN THE PEGASUS DATA, POSITION BY POSITION (3 

C @ (3 

C 

C FIRST ZERO OUT THE P MATRIX 
C 

, DO 20 I = 1,NT 
DO 20 J = 1,NT 
P(I,J) = 0. ODO 
20 CONTINUE 
N = 0 

25 READd,*, END = 45) PRESS, ZP, (TRANS( I) , I=1,NT) ,SV(N+1) 

NNT = N*NT 
DO 30 I = 1, NT 
LB(NOT +!,!)= TRANS(I) 

30 CONTINUE 

CALCULATE THE APPROXIMATE POSITION OF PEGASUS 
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CALL APPR0(TRANS,X,Y,Z,SV(N+1) ,NT.AUX,XP,YP,ZP) 

C 

X0(N*3 + 1) = XP 

X0(N*3 + 2) = YP 

X0(N*3 + 3) = ZP 

WRITE(2 32) PRESS,ZP,(TRANS(I) ,1=1 ,NT) ,SV(N+1) ,XP,YP,ZP 

32 FORMAT(' ' , 3X,F7. 2 , 2X ,F7. 2 ,5X ,F7. A , 3( 6X ,F7. 4) , 13X ,F8. 3 , 7X , 

1 3(2X,F7. 1)) 

DO ERROR PROPAGATION AND FORM THE P MATRIX FOR THE FIRST PEGASUS 
POSITION. STORE THIS IN MATRIX BIGP. THE ERROR PROPAGATION 
SUBROUTINE MUST BE INSERTED HERE AND THE NEXT FEW LINES OF CODE 
MODIFIED ACCORDINGLY. 

Ex.ExSS-Ex.Ex.&.&.&.StExEt&.&.&.Ex.&-&S.S.Ex.&.S.Ex&.SLS.Ex.&-&-St&.StSt&.SLS.Ex&.^Ex.&. 

FOR INITIAL PROGRAM TESTING, INSERT DIAGONAL VALUES ONLY 

DO 35 I = 1, NT 
P(I,I) = 1. ODO/0. 0005D0**2 
35 CONTINUE 

DO 40 I = 1,NT 
DO 40 J = 1,NT 
BIGP(I,NNT + J) = P(I,J) 

40 CONTINUE 
N = N + 1 
GO TO 25 

45 CONTINUE 

CALCULATE PROGRAM CONSTANTS 
NMAX = N 

NP = N*3 + (NT’'^3) 

NOBS = N*NT 
NV = (6*NT) * N 
NN = NP*(NP+l)/2 
NT3= NT* 3 
ITER = 0 

WRITE(2,46) NMAX, NOBS, NP 

46 FORMATC ' 1' ,///,5X, 'ACTUAL // OF PEGASUS POSITIONS =’,I5,//, 

1 5X, 'TOTAL it OF OBS. ON PEGASUS =',I5,//, 

2 5X, 'TOTAL it OF PARAMETERS =',I5) 

C @ 

C (§ SET UP THE MATRICES NEEDED FOR THE ADJUSTMENT, BEGINNING 

C @ WITH THOSE ASSOCIATED WITH THE "OBSERVED" PARAMETERS. 

C @ MOST OF THESE NEED TO BE ZEROED OUT AT THIS POINT. 

C @ 

C @ THEN SET UP THE NORMAL MATRIX BLOCK BY BLOCK, EACH BLOCK 

C @ CORRESPONDING' TO A NEW PEGASUS POSITION 

C @ 
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ZERO OUT MATRICES 

DO 50 I = 1,NP 
LB2(I) = 0. ODO 
P2( I) = 0. ODO 
L2(I) = O.ODO 
XHAT(I) = 0. ODO 
50 CONTINUE 

INSERT THE NON-ZERO ELEMENTS 
J = N*3 

DO 55 I = 3,J,3 
P2(I) = 1. 0D0/SDZP**2 
LB2(I)= XO(I) 

55 CONTINUE 
C 



J = N*3 


4- 


1 


K = 1 






DO 58 I 


= 


J,NP,3 


XO(I) 


= 


X(K) 


P2(I) 


= 


1. 0D0/SDX(K)**2 


LB2(I) 


= 


XO(I) 


X0(I+1) 


= 


Y(K) 


P2(I+1) 


= 


1. 0D0/SDY(K)**2 


LB2(I+1) 


= 


X0(I+1) 


X0(I+2) 


= 


Z(K) 


P2(I+2) 


= 


1. 0D0/SDZ(K)**2 


LB2(I+2) 


= 


X0(I+2) 


K = K + 


1 





58 CONTINUE 

WRITE(2,800)(LB(I,1) ,I=l,NOBS) 

800 FORMAT( ' l' ,///,50X,'LB MATRIX' ,// ,8( IX, 10(F10. 7 ,2X)/) , 1X.2F10. 7) 
WR1TE(2,801)(X0(I),I=1,NP) 

801 F0RMAT('0' ,///,50X,'X0 VECTOR* , // , 8( IX, 10(F10. 4 , 2X)/) , IX ,F10. 4) 
WRITE(2,802)(P2(I) ,I=1,NP) 

802 FORMATC ' O’ ,/// ,50X, ’ P2 VECTOR* ,//,8( IX , 10(F10. 4 ,2X)/) , IX ,F10. 4) 
WRITE(2,803)(LB2( I) ,I = 1,NP) 

803 FORMAT(*0* , /// , 50X , * LB2 VECTOR* , // , 8( IX , 10(F10. 4, 2X)/) , 1X,F10. 4) 
C 

C (3 (3 

C @ THE ITERATIVE PROCESS BEGINS FROM HERE. BEGIN BY ZEROING (3 

C (3 OUT THE NORMAL MATRIX AND THE TWO COMMON BLOCK MATRICES @ 

C (3 0 

C 

60 DO 62 I = 1,NT3 
SPL(I,1) = 0. ODO 
DO 62 J = 1,NT3 
SPS(I,J) = 0. ODO 
62 CONTINUE 

DO 65 I = 1,NN ■ 

ANORM( I) = 0. ODO 
65 CONTINUE 
C 
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K = 1 

J = N*3 + 1 
DO 66 I = J,NP,3 
X(K) = XO(I) 
Y(K) = X0(I+1) 
Z(K) = X0(I+2) 
K = K+1 
66 CONTINUE 



BEGIN THE FORMATION OF THE A AND LI MATRICES 

DO 100 K = l.NMAX 
JJ = (K-1)*NT 
DO 68 I = 1,NT 
TRANS(I) = LB(JJ + 1,1) 

DO 68 J = 1,NT 

P(I,J) = BIGP(I,JJ + J) 

68 CONTINUE 

KK = (K-l)*3 
XP = X0(KK + 1) 

YP = X0(KK + 2) 

ZP = X0(KK + 3) 



DO 70 I = 1,NT 

L1(JJ + 1,1) = (DSQRT((XP-X(I))**2 + (YP-Y(I))**2 + (ZP-Z( I) )**2) 
1 /SV(K)) - TRANS(I) 

L(I,1) = L1(JJ + 1,1) 

70 CONTINUE 

ZERO OUT THE A AND S SUBMATRICES 

DO 80 I = 1,NT 
DO 75 J = 1,3 
A(I,J) = 0. ODO 
75 CONTINUE 

DO 80 J = 1,NT3 
S(I,J) = O.ODO 
80 CONTINUE 

CALL SETUP(X,Y,Z,XP,YP,ZP,SV(K) ,K ,NT,NV ,NT3 ,NMAX , A, S , ICOL, IROW , 

1 VALUE) 

CALL NORM(A,S,P,L,K,NT,NT3,NN,NP,NMAX,ATP,ATPA,ATPS,ATPL,STP, 

1 STPS,STPL,ANORM,XHAT) 

100 CONTINUE 

WRITE(2,804)(L1(I,1) ,I=l,NOBS) 

804 FORMATCi’ ,///,50X,'L1 MATRIX' , // ,8( IX, 10(F10. 7 ,2X)/) , 1X,2F10. 7) 

PLACE SPS AND SPL INTO THE NORMAL MATRIX AND THE COLUMN VECTOR 
RESPECTIVELY. 

DO no I = 1,NT3 
NCOL = (NP + 1) - I 
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11 = NCOL*(NCOL + l)/2 + I 
DO 110 J = I,NT3 

12 = II - J 

ANORM(I2) = SPS(NT3-J+1,NT3-I+1) 
no CONTINUE 

13 = NMAX*3 + 1 
DO 120 I = 13, NP 
XHAT(I) = SPL(I-I3+1,1) 

120 CONTINUE 

SET UP THE L2 MATRIX AND ADJUST THE DIAGONAL ELEMENTS OF THE 
NORMAL MATRIX FOR THE INFLUENCE OF P2 

CALL P2L2 ( P2 , LB 2 , XO , ANORM , L2 , XHAT , NP , NN ) 

(a (3 

(a SOLVE THE NORMAL EQUATIONS AND COMPUTE THE RESIDUALS @ 

(a @ 

^(3(§(3(3(300(9(3(9@(3(9(9(3@00@(^@(^(3(3(300@(50(3(3(2(3@(3(3@0(9(9(3^@(^00(S(3@@@l(ci(3(3(3@(3 

call a user supplied routine tu solve the normal equations. 

(SEE APPENDIX A, SECTION 0.) 

DO 135 I = 1,NP 
XHAT(I) = -XHAT(I) 

XA(I) = XO(I) + XHAT(I) 

XO(I) = XA(I) 

135 CONTINUE 

ITER = ITER + 1 
WRITE(2,140) 

140 FORMAT! 'O' ,///,5X, 'ADJUSTED PEGASUS POSITIONS (WITH THE ', 

1 'CORRECTIONS TO THE APPROX. POSITIONS IN PARENTHESES)',//, 18X, 

2 'X' ,18X,'Y' ,18X,'Z’ ) 

J = NP - NT3 

DO 145 1=1, J, 3 

WRITE(2,142) XA(I) ,XHAT(I) ,XA( I + l) ,XHAT( I + l) ,XA(I+2) , 

1 XHAT(I+2) 

142 FORMAT( ' ' ,10X,F9. 3, 2X , ' ( ' , F6. 1 , ' ) ‘ , 2( 2X , F9. 3 , 2X, ' ( ’ ,F6. 1, ' )' )) 

145 CONTINUE 
J = J + 1 
WRITE(2,146) 

146 FORMATC ’ ,//,5X,' ADJUSTED TRANSPONDER POSITIONS (WITH’, 

1 'CORRECTIONS TO THE APPROX. POSITIONS IN PARENTHESES)') 

DO 148 I=J,NP,3 

WRITE(2,142) XA(I) ,XHAT(I) ,XA(I+1) ,XHAT( I+l) ,XA( 1+2) , 

1 XHAT(I+2) 

148 CONTINUE 

TEST THE PARAMETER CORRECTIONS TO SEE IF ITERATION IS REQUIRED 

I FLAG = 0 
DO 150 I =1,NP 

IF(DABS(XHAT(I)).GT.TOL)IFLAG = IFLAG + 1 
150 CONTINUE 
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C COMMENT OUT THE NEXT THREE STATEMENTS DURING PROGRAM DEVELOPMENT 
IF( IFLAG. EQ. 0) GO TO 155 

IF( ITER. LT. NITER) GO TO 60 
C 

COMPUTE THE RESIDUALS 

155 CALL RESID(VALUE,IC0L,IR0W,BIGP,L1,XHAT,NMAX,NT,NV,NP,N0BS, 

1 P.VTP.V.SIG.Vl ,SIG0) 

WRITE(2 160) 

160 FORMATC ' 1 ',///, lOX, 'OBSERVATIONAL RESIDUALS AFTER THE ADJUSTMENT', 
1 //,5X, 'TRANSPONDER 1 TRANSPONDER 2 TRANSPONDER 3' ,4X, 

2 'TRANSPONDER 4') 

DO 170 I = 1, NOBS, NT 
WRITE(2, 165)(V1( I+K-1) , K=1,NT) 

165 FORMATC^ ' , 9X , F8. 4 , 3( 9X , F8. 4 ) ) 

170 CONTINUE 

WRITE(2 175) SIGO 

175 FORMATC* ' ,/// ,5X, 'THE A POSTERIORI VARIANCE OF UNIT WEIGHT =’ , 

1 F8.4) 

*THE FOLLOWING WRITE STATEMENT WAS ADDED BY M. HASKELL TO LOOK 
*AT THE VARIANCE-COVARIANCE MATRIX FOR POSITION COORDINATES. 

*THE CALL STATEMENT FOR SUBROUTINE MAT WAS ADDED TO ACCESS ONLY 
*THE INFORMATION NEEDED TO CALCUL.4TE THE VARIANCES AND COVARIANCES 
*0F THE VELOCITIES. 

WRITE(2,*)(AN0RM(II) ,II=1,NN) ^PROVIDES (ATPA)"^ IN VECTOR FORM 
CALL MAT(AN0RM,NMAX) 



(a (a 
@ COMPUTE VARIOUS STATISTICAL QUANTITIES INCLUDING THE @ 
(§ VELOCITY OF PEGASUS. (? 
0 0 






STOP 

END 



0 0 

0 SUBROUTINES USED IN THE PROGRAM 0 

0 0 

(a(a(a^@@0(a(^0^(^(a^^(a(^(a(a(^(a{a^(a@@(^(^(a{a(a@@0(a(a(^(a 



SUBROUTINE APPRO( TRANS , X , Y , Z , SV , NT , AUX , XP , YP , ZP) 

IMPLICIT REAL*8(A-H,0-Z) 

THIS SUBROUTINE COMPUTES THE APPROXIMATE COORDINATES FOR PEGASUS 
INPUT; TRANS - VECTOR OF TRANSPONDER TIME DELAYS FOR THE DESIRED 
PEGASUS POSITION. 
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X,Y,Z - ARRAYS CONTAINING THE COORDINATES OF THE TRANSPONDERS 
SV - VELOCITY OF SOUND IN WATER 
NT - NUMBER OF TRANSPONDERS 
OUTPUT: XP,YP,ZP - COORDINATES OF PEGASUS 

DIMENSION TRANS(NT), X(NT) , Y(NT) , Z( NT) , AUX( 2) 

DO 10 I = 1,NT 
TRANS(I) = TRANS(I)*SV 
10 CONTINUE 

DO 20 I = 1,2 

AUX(I) = (TRANS(I + 1)**2 - ( ZP-Z( 1+1) )**2 - X(I + 1)**2 - Y(I + 1)**2) 

1 - (TRANS(I)**2 - (ZP-Z(I))**2 - X(I)**2 - Y(I)**2) 

20 CONTINUE 

YP = 0. 5D0*(AUX(1)/(X(1)-X(2)) - AUX( 2)/(X( 2) -X( 3) ) )/( ( Y( 1) -Y( 2) ) 

1 /(X(l)-X(2)) - (Y(2)-Y(3))/(X(2)-X(3))) 

XP = 0. 5D0*(AUX( 1)/(Y(1)-Y(2)) - AUX( 2)/( Y( 2) -Y( 3) ) )/( (X( 1) -X( 2) ) 

1 /(Y(l)-Y(2)) - (X(2)-X(3))/(Y(2)-Y(3))) 

DO 30 I = l.NT 
TRANS(I) = TRANS(I)/SV 
30 CONTINUE 
RETURN 
END 

SUBROUTINE SETUP(X , Y , Z ,XP , YP , ZP, SV,N ,NT,NV,NT3 ,NMAX, A,S, ICOL, IROW, 
1 VALUE) 

IMPLICIT REAL*8(A-H,0-Z) 

THIS MATRIX SETS UP THE A AND S MATRICES NEEDED FOR THE BLOCK 
BY BLOCK FORMATION OF THE NORMAL MATRIX. 

INPUT: X(I), I = l.NT X COORDINATES FOR THE TRANSPONDERS 

Y(I), I = 1,NT Y 
Z(I), I = 1,NT Z 

XP.YP.ZP, APPROXIMATE COORDINATES FOR PEGASUS. 

SV = VELOCITY OF SOUND THROUGH WATER. 

N = NTH POSITION OF PEGASUS. 

NT, NV, NT3 SAME AS IN THE MAIN PROGRAM. 

NMAX = ACTUAL if OF PEGASUS POSITIONS IN THIS DROP 
OUTPUT: A = THE PORTION OF THE A MATRIX CORRESPONDING TO THE 
NTH PEGASUS POSITION. 

S = AS ABOVE, BUT THE OUTER BAND OF THE A MATRIX. 

VALUE = THE NON ZERO if IN THE A MATRIX CORRESPONDING 

TO THE COLUMN if HELD IN ICOL AND THE ROW if HELD 
IN IROW. 

DIMENSION X(NT), Y(NT), Z(NT), A(NT,3), S(NT,NT3), VALUE(NV) 
INTEGER*2 ICOL(NV), IROW(NV) 

NROW = (N-1)*NT 
NCOL = (N-l)*3 
NCOLl = NMAX*3 
ICOUNT = (N-1)*6*NT 
DO 20 I = l.NT • 

J = 1 

K = (I-l)*3 + 1 

DIST = DSQRT((XP - X(I))**2 + (YP - Y(I))**2 + (ZP - Z(I))**2) 



C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 



C 
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A(I,J) = (XP - X(I))/(DIST*SV) 

I COUNT = I COUNT + 1 
VALUE ( I COUNT) = A(I,J) 

IROW(ICOUNT) = NROW + I 
ICOL(ICOUNT) = NCOL + J 
A(I,J+1) = (YP - Y(I))/(DIS'P'^SV) 

ICOUNT = ICOUNT + 1 
VALUE(ICOUNT) = A(I,J+1) 

IROW( ICOUNT) = NROW + I 
ICOL( ICOUNT) = NCOL + J+1 
A(I,J+2) = (ZP - Z(I))/(DIS'pvsv) 

ICOUNT = ICOUNT + 1 
VALUE( ICOUNT) = A(I,J+2) 

IROW( ICOUNT) = NROW + I 
ICOL( ICOUNT) = NCOL + J+2 
S(I,K) = -A(I,J) 

ICOUNT = ICOUNT + 1 
VALUE( ICOUNT) = S(I,K) 

IROW( ICOUNT) = NROW + I 
ICOL( ICOUNT) = NCOLl + K 
S(I,K+1) = -A(I,J+1) 

ICOUNT = ICOUNT + 1 
VALUE( ICOUNT) = S(I,K+1) 

I ROW( ICOUNT) = NROW + I 
ICOL( ICOUNT) = NCOLl + K+1 
S(I,K+2) = -A(I,J+2) 

ICOUNT = ICOUNT + 1 
VALUE( ICOUNT) = S(I,K+2) 

IROW( ICOUNT) = NROW + I 
ICOL( ICOUNT) = NCOLl + K+2 
20 CONTINUE 
WRITE(2,100) 

100 FORNATC ' 1' ,///,15X,'A MATRIX AND THEN THE S MATRIX, BLOCK BY', 

1 ’ BLOCK') 

WRITE(2,101)((A(I,J) ,J=1,3) ,1 = 1 ,4) 

101 FORMATC ', T2,4(2X,3(F10. 8,3X) ,/)) 

WRITE(2,102)((S(I,J) ,J=1, 12) ,1=1,4) 

102 FORMATC T2 ,4( 12(F10. 7 , IX) , / ) ) 

RETURN 

END 

SUBROUTINE NORM( A,S , P , L,N ,NT,NT3 ,NN,NP ,NMAX , ATP ,ATPA, ATPS ,ATPL, 
1 STP,STPS,STPL,ANORM,XHAT) 

IMPLICIT REAL*8(A-H,0-Z) 

THIS SUBROUTINE FORMS THE NORMAL EQUATIONS, BLOCK BY BLOCK 
INPUT: A MATRIX FROM SUBROUTINE SETUP 

g n It It tt I 

P THE WEIGHT MATRIX CORRESPONDING TO THIS BLOCK 
L THE L MATRIX 

N, NT, NT3, AS FOR SUBROUTINE SETUP 
NN, NP, NMAX, AS FOR MAIN PROGRAM 
OUTPUT: A NUMBER OF AUXILLARY MATRICES USED IN THE MATRIX 

MULTIPLICATIONS, IE. , ATP, ATPA, ATPS, ATPL,STP,STPS, 
AND STPL. 

MATRICES SPS AND SPL WHICH COME IN VIA THE COMMON 
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BLOCK ARE UPDATED. 

ANORM - THE NORMAL MATRIX (IN ARRAY FORM) 
XHAT - THE ATPL VECTOR 



DIMENSION A(NT,3), S(NT,NT3) ,P(NT,NT) ,ATP(NT,3) , ATPA(3,3), 

1 ATPS(3,NT3), ATPL(3,1), STP(NT3,NT), STPS(NT3,NT3) ,STPL(NT3, 1) , 

2 ANORM(NN) ,XHAT(NP) 

DOUBLE PRECISION L(NT,1) 

COMMON SPS(12,12), SPL(12,1) 

PERFORM THE MATRIX MULTIPLICATIONS 

CALL ATB(A,P,ATP,NT,3,NT) 

CALL AB(ATP,A,ATPA,3,NT,3) 

CALL AB(ATP,S,ATPS,3,NT,NT3) 

CALL AB(ATP,L,ATPL,3,NT,1) 

CALL ATB(S,P,STP,NT,NT3,NT) 

CALL AB(STP,S,STPS,NT3,NT,NT3) 

CALL AB(STP,L,STPL,NT3,NT, 1) 

WRITE(2 100)((ATPA(I,J),J=1,3),ATPL(I,1),I=1,3) 

100 FORMATC' ',/,15X, 'ATPA BLOCK* , lOX ATPL VECTOR*,//, 

1 3(2X,3(D10. 3,2X) , 4X,D10.3,/)) 
WRITE(2,110)((STPS(I,J),J=1,12),I=1,12) 

110 FORMATC* *,/,15X, ' STPS BLOCK* ,// , 

1 T2,12(12(D10. 3,1X),/)) 

PLACE ATPA AND ATPL IN THEIR APPROPRIATE POSITION IN EITHER THE 
NORMAL MATRIX OR THE COLUMN VECTOR 

NCOL = (N-l)*3 + 1 

II = NCOL*(NCOL + l)/2 
ANORM(Il) = ATPA(1,1) 

III = II + NCOL 
ANORM(Ill) = ATPA(1,2) 

ANORMdll + 1) = ATPA(2,2) 
nil = 111 + 1 + NCOL 
ANORM(Illl) = ATPA(1,3) 

ANORMdlLl + 1) = ATPA(2,3) 

AN0RM(I111 + 2) = ATPA(3,3) 

DO 10 I = 1,NT3 
INT = (NMAX*3) + I 
J = INT’KINT - l)/2 + 1 + (N-l)*3 
ANORM(J) = ATPS(1,I) 

AN0RM(J+1) = ATPS(2,I) 

ANORM(J+2) = ATPS(3,I) 

10 CONTINUE 



DO 20 I = 1,3 
XHAT(NCOL+I-l) =ATPL(I,1) 
20 CONTINUE 



UPDATE SPS AND SPL 



DO 30 I = 1,NT3 
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SPL(I.l) = SPL(I,1) + STPL(I.l) 
DO 30 J = 1.NT3 

SPS(I.J) = SPS(I,J) + STPS(I,J) 
30 CONTINUE 
RETURN 
END 



SUBROUTINE ATB( A , B . R . L.M ,N) 

IMPLICIT REAL’'-8(A-H,0-2) 

THIS SUBROUTINE COMPUTES THE MATRIX PRODUCT A’B 
INPUT; MATRIX A (L X M) 

MATRIX B (L X N) 

OUTPUT MATRIX R (M X N) 

DIMENSION A(L,M) ,B(L,N) ,R(M,N) 

DO 10 I = 1,M 
DO 10 J = 1,N 
R( I ,J) = 0. ODO 
DO 5 K = 1 , L 

R(I ,J) = R(I ,J) + A(K,I)*B(K,J) 

5 CONTINUE 
10 CONTINUE 
RETURN 
END 



SUBROUTINE AB( A , B , R , L,M , N) 

IMPLICIT REAL*8(A-H,0-Z) 

FORM THE MATRIX PRODUCT R = AB. 

THE MATRICES A AND B ARE RETURNED UNCHANGED 

DIMENSION A(L,M) ,B(M,N) ,R(L,N) 

DO 5 I = I,L 
DO 5 J = 1,N 
R(I,J) = 0. ODO 
DO 5 K = l.M 

R(I,J) = R(I,J) + A(I,K)*B(K,J) 

5 CONTINUE 
RETURN 
END 



SUBROUTINE RES I D( VALUE , I CO L , I ROW , B I GP , L 1 , XHAT , NMAX , NT , NV , NP , NOB S , 
1 P,VTP,V,SIG,V1,SIG0) 

IMPLICIT REAL*8(A-H,0-Z) 

COMPUTES THE RESIDUALS ON ALL THE TRANSPONDER TIME DELAY 
OBSERVATIONS AND THE A POSTERIORI VARIANCE OF UNIT WEIGHT. 

INPUT; VALUE, ICOL.IROW,- MATRICES DERIVED IN SUBROUTINE SETUP 
BIGP - THE FULL WEIGHT MATRIX FOR ALL OBSERVATIONS 
LI - THE VECTOR OF "COMPUTED-OBSERVED" OBSERVATIONS 
XH.AT - THE LEAST SQUARES SOLUTION VECTOR 
NMAX, NT, NV.NP, NOBS - AS IN THE MAIN PROGRAM 
AUXILLARY MATRICES : P, VTP , V,SIG 
OUTPUT; VI - THE VECTOR OF RESIDUALS 

SIGO - THE A POSTERIORI VARIANCE OF UNIT WEIGHT 
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DIMENSION VALUE(NV) ,XHAT(NP) , Vl(NOBS) ,BIGP(NT,NOBS) ,P(NT,NT) , 
I VTP(1,NT),SIG(1,1),V(NT,1) 

DOUBLE PRECISION LI(NOBS) 

INTEGER*2 ICOL(NV), IROW(NV) 

DO 10 I = l.NOBS 
V1(I) = O.ODO 
10 CONTINUE 
IC = 1 

15 DO 20 I = 1,NV 

IF(IROW(I).NE. IC) GO TO 18 

16 Vl(IC) = Vl(IC) + VALUE(I)*XHAT(ICOL(I)) 

GO TO 20 

18 Vl(IC) = Vl(IC) + Ll(IC) 

IC = IC + 1 

IF(IC. LE. NOBS) GO TO 16 
20 CONTINUE 

Vl(NOBS) = Vl(NOBS) + Ll(NOBS) 

SIGO = 0. ODO 
1C = 0 

DO 40 I = l.NMAX 
DO 30 J = 1,NT 
V(J,1) = V1(IC*NT + J) 

DO 30 K = 1,NT 

P(J,K) = BIGP(J,IC*NT + K) 

30 CONTINUE 

CALL ATB(V,P,VTP,NT,1,NT) 

CALL AB(VTP,V,SIG,1,NT,1) 

SIGO = SIGO + SIG(1,1) 

IC = IC + 1 
40 CONTINUE 

SIGO = SIGO/FLOAT(NOBS-NMAX*3) 

RETURN 

END 



SUBROUTINE P2L2( P2 , LB2 , XO , ANORM , L2 , XHAT , NP , NN ) 

IMPLICIT REAL*8(A-H,0-Z) 

THIS SUBROUTINE COMPUTES THE L2 MATRIX AND ADDS IT TO THE COLUMN 
VECTOR. IN ADDITION IT INCREMENTS THE DIAGONAL ELEMENTS OF THE 
NORMAL MATRIX TO INCLUDE THE INFLUENCE OF P2 

DIMENSION P2(NP), XO(NP) .ANORM(NN) ,XHAT(NP, 1) 

DOUBLE PRECISION L2(NP), LB2(NP) 

DO 10 I = 1,NP 
L2(I) = XO(I) - LB2(I) 

INT = I*(I+l)/2 

ANORM(INT) = ANORM(INT) + P2(I) 

P2L = P2(1)*L2(I) 

XHAT(I.I) = XHAT(I,1) + P2L 
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10 CONTINUE 
RETURN 
END 



C 

C 

C 

C 

C 



C 



C 



C 



C 



C 



C 



C 



C 



SUBROUTINE MAT(C,NMAX) 

*NOT PART OF J HANNAH'S ORIGINAL PROGRAM. ADDED BY M. HASKELL. 

THIS PROGRAM READS APPROPRIATE VALUES FROM THE VARIANCE-COVARIANCE 
MATRIX FILE (ANORM) AND PLACES THEM IN THE MATRICES NEEDED TO COM- 
PUTE VARIANCES OF THE CURRENT VELOCITIES. 

IMPLICIT REAL*8 (A-H.O-Z) 

DOUBLE PRECISION C(l) 

DIMENSION C(4A4153) 

I(K) = K*(K+l)/2 
K = 1 
ROW = K 

WRITE (2,*)C(I(K)), C(I(K+1)-1), C(I(K+2)-2) 

WRITE (2,*)C(I(K+1)-1) , C(I(K+1)), C(l(K+2)-l) 

WRITE (2,*)C(I(K+2)-2) , C(I(K+2)-l), C(I(K+2)) 

DO 15 J = 1,NMAX+1 

WRITE (2,*)C(I(K+3)-3), C(I(K+4)-4), C(I(K+5)-5) 

WRITE (2,*)C(I(K+3)-2) , C(I(K+4)-3), C(I(K+5)-4) 

WRITE (2,*')C(I(K+3)-l) , C(I(K+4)-2), C(I(K+5)-3) 

WRITE (2,*)C(I(K+3)) , C(I(K+4)-l), C(I(K+5)-2) 

WRITE (2,*)C(I(K+4)-l), C(I(K+4)), C(I(K+5)-l) 

WRITE (2,*)C(I(K+5)-2) , C(I(K+5)-l), C(I(K+5)) 

K = K+3 
15 CONTINUE 



WRITE (2,*)C(I(K+3)-6) , C(I(K+4)-7), C(I(K+5)-8) 
WRITE (2,*)C( I(K+3)-5) , C(I(K+4)-6), C(I(K+5)-7) 
WRITE (2,*)C(I(K+3)-4) , C(I(K+4)-5), C(I(K+5)-6) 

WRITE (2,*)C(I(K+3)-3) , C(I(K+4)-4), C(I(K+5)-5) 
WRITE (2,*)C(I(K+3)-2), C(I(K+4)-3), C(I(K+5)-4) 
WRITE (2,*)C( I(K+3)-l) , C(I(K+4)-2), C(I(K+5)-3) 

WRITE (2,*)C(I(K+3)) , C(I(K+4)-l), C(I(K+5)-2) 
WRITE (2,*)C( I(K+4)-l) , C(I(K+4)), C(I(K+5)-l) 
WRITE (2,*)C( I(K+5)-2) , C(I(K+5)-l), C(I(K+5)) 



K = K+3 

WRITE (2,*)C(I(K+3)-9) , 
WRITE (2,*)C(I(K+3)-8) , 
WRITE (2,*)C(I(K+3)-7), 



C(I(K+4)-10) , C(I(K+5)-ll) 
C(I(K+4)-9), C(I(K+5)-10) 
C(I(K+4)-8), C(I(K+5)-9) 



WRITE (2,*)C(I(K+3)-6) , C(I(K+4)-7), 
WRITE (2,*)C( I(K+3)-5) , C(I(K+4)-6), 
WRITE (2,*)C(I(K+3)-4), C(I(K+4)-5), 



C(I(K+5)-8) 
C(I(K+5)-7) 
C( I(K+5)-6) 



WRITE (2,*)C(I(K+3)-3), C(I(K+4)-4), 
WRITE (2,*)C(I(K+3)-2), C(I(K+4)-3), 
WRITE (2,*)C(I(K+3)-l), C(I(K+4)-2), 



C(I(K+5)-5) 
C( I(K+5)-4) 
C(I(K+5)-3) 
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WRITE (2,’'')C(I(K+3)) , C( I ( K-+ /. ) * 1 ) . C(I(K + 5'l-2) 
WRITE (2 ,>v)C(I(K+ 4)-1) , C(I(K+4)), C(I(,K+5)-l) 
WRITE (2,*)C(I(K+5)-2) , C(I(K+5)-l), C(I(K+5)) 

K = 36 

DO 15 J = 1,6 

WRITE (2,*)C(I(K+3)-3) , C(I(K+4)-4), C(I(K+5)-5) 
WRITE (2,*)C(I(K+3)-2) , C(I(K+4)-3), C(I(K+5)-4) 
WRITE (2,*)C( I(K+3)-l) , C(I(K+4)-2), C(I(K+5)-3) 

WRITE (2,*)C(I(K+3)) , C(I(K+4)-l), C(I(K+5)-2) 
WRITE (2,*)C(I(K+4)-l) , C(I(K+4)), C(I(K+5)-l) 
WRITE (2,*)C(I(K+5)-2) , C(I(K+5)-l), C(I(K+5)) 

K = K+36 
15 CONTINUE 
RETURN 
END 
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CURRENT VELOCITY VARIANCE-COVARIANCE PROGRAM 



This program computes the variance-covariance matrix of 
the current velocities by a procedure described in Chapter VI, 
Section D. It uses data from the Pegasus position variance- 
covariance matrix produced by the Fortran program Pegasus as 
described in Appendix A. 
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C THIS PROGRAM USES DATA FROM THE VARIANCE-COVARIANCE MATRIX PRODUCED 
C BY THE PEGASUS FORTRAN PROGRAM. IT COMPUTES THE VARIANCE-COVARI - 
C ANCE MATRIX OF THE CURRENT VELOCITIES. 

C 

IMPLICIT REAL*8(A-H,0-Z) 

REAL SP1(3,3), SP12(3,3), SP2( 3 , 3) , V( 3 , 3) ,SU , SV.SW 
C 

N = 1 

DO 100 I = 1,3 

READ(1,*)(SP1(I,J),J = 1,3) 

C WR1TE(2,120)(SP1(I,J),J = 1,3) 

120 FORMATC' ' , 3X ,F10. 4 , 3X,F10. 4 , 3X ,F10. A) 

100 CONTINUE 
C 

10 DO 200 I = 1,3 

READ(1,*,END=45)(SP12(I,J),J = 1,3) 

C WRITE(2,120)(SP12(I,J),J = 1,3) 

200 CONTINUE 

DO 300 I = 1,3 

READ( 1,*)(SP2( I ,J) ,J = 1,3) 

C WR1TE(2,120)(SP2(I,J),J = 1,3) 

300 CONTINUE 
C 

CC WRITE(2,15)N 

15 FORMAT(' '/,5X,' VELOCITY VARIANCE-COVARIANCE MATRIX FOR ', 

I'POSITION ,13,/) 

DO 400 I = 1,3 
DO 350 J = 1,3 

V(I,J)=((SP1(I,J)+SP2(I,J)-SP12(I,J)-SP12(J,I))/(16**2) 

350 CONTINUE 

C WRITE(2.120)(V(I,J),J = 1,3) 

400 CONTINUE 

SU = SQRT(V(1,1)) 

SV = SQRT(V(2,2)) 

SW = SQRT(V(3,3)) 

CC WRITE(2,130)SU, SV, SW 

CC130 FORMATC ' ',/ ,3X, ’ SIGMA U =' ,F10. 4 ,3X, ' SIGMA V =',F10. 4,3X, 

CC 1 'SIGMA W =' ,F10. 4/) 

WRITE(2,600)N,SU,SV,SW 

600 FORMATC ' ',10X,'P0S ' , 13 , 3X,F7. 4 , 3X,F7. 4 , 3X ,F7. 4) 

DO 450 I = 1,3 
, DO 460 J = 1,3 
SP1(1,J) = SP2(I,J) 

460 CONTINUE 

WRITE(2,140)(SP1(I,J),J=1,3) 

140 FORMATC* ' , 3X,F10. 4 , 3X ,F10. 4 , 3X ,F10. 4) 

450 CONTINUE 
N = N +1 
GO TO 10 
45 CONTINUE 
STOP 
END 
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